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Abstract 

A  recently  developed  class  of  models  incorporating  the  cyton  model  of  population  generation  structure  into 
a  conservation-based  model  of  intracellular  label  dynamics  is  reviewed.  Statistical  aspects  of  the  data  collection 
process  are  quantified  and  incorporated  into  a  parameter  estimation  scheme.  This  scheme  is  then  applied  to 
experimental  data  for  PHA-stimulated  CD4+  T  and  CD8+  T  cells  collected  from  two  healthy  donors.  This 
novel  mathematical  and  statistical  framework  is  shown  to  form  the  basis  for  accurate,  meaningful  analysis  of 
cellular  behavior  for  a  population  of  cells  labeled  with  the  dye  CFSE  and  stimulated  to  divide. 


1  Introduction 

Dating  at  least  as  far  back  as  the  work  of  Bell  and  Anderson  [14]  in  1967,  mathematical  models  have  been 
proposed  which  attempt  to  describe  the  biophysical  processes  involved  in  cell  division.  These  models  range  in 
scope  and  scale  from  phenomenological  descriptions  of  population  growth  to  mechanistic  models  of  subcellular 
mechanics.  One  particularly  important  class  of  mathematical  models  is  that  in  which  the  behaviors  of  individual 
cells  are  linked  in  a  meaningful  way  to  population  level  characteristics.  This  class  of  models  has  applications  in 
quantitative  descriptions  of  the  immune  system,  where  the  behavior  of  individual  cells  can  vary  widely  across  the 
population  but  in  which  the  immune  response  (understood  to  be  the  net  result  of  actions  of  all  relevant  cells  in 
the  system)  is  much  more  predictable  [43].  In  fact,  a  quantitative  description  of  the  ‘cellular  calculus’  [25]  by 
which  cells  send,  receive,  and  respond  to  intra-  and  extracellular  stimuli  is  in  many  respects  an  open  problem  in 
immunology. 

In  the  past,  the  major  limitation  of  mathematical  models  linking  cellular  and  population- level  behavior  has 
been  the  difficulty  in  obtaining  data  with  which  to  validate  the  models.  Generally,  it  is  possible  to  observe  the 
behavior  of  a  small  number  of  single  cells  quite  carefully  in  isolation,  or  to  observe  a  large  number  of  cells 
in  aggregate.  More  recently,  the  intracellular  dye  carboxyfluorescein  succinimidyl  ester  (CFSE)  [37]  for  use  in 
flow-cytometric  proliferation  assays  in  vitro  has  emerged  as  a  powerful  experimental  technique  for  the  study  of 
dividing  cells.  Because  the  dye  emits  bright,  approximately  uniform,  and  long-lasting  fluorescent  labeling  of  a 
population  of  cells  and  is  approximately  evenly  partitioned  during  cell  division,  the  dye  provides  a  useful  surrogate 
for  the  number  of  divisions  a  cell  has  undergone.  Individual  cells  can  be  assessed  by  a  flow  cytometer,  which 
can  simultaneously  measure  additional  properties  of  cells  such  as  size,  internal  complexity,  cell  surface  marker 
expression,  levels  of  cytokine  secretion,  etc. 

When  a  population  of  cells  is  measured,  the  individual  CFSE  fluorescence  intensity  measurements  can  be 
placed  into  a  histogram  as  in  Figures  1  and  2.  Each  ‘peak’  in  the  histogram  represents  a  cohort  of  cells  having 
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completed  the  same  number  of  divisions.  When  such  measurements  are  made  sequentially  in  time,  one  obtains 
information  on  the  dynamic  response  of  the  population  of  cells  to  a  stimulus.  As  such,  CFSE-based  flow  cytometric 
analysis  is  a  promising  tool  for  the  study  of  cell  division  and  division-linked  changes.  The  ultimate  goal  for  the 
quantitative  analysis  of  CFSE  data  (in  particular,  as  it  relates  to  studies  of  the  immune  system)  is  to  incorporate 
fundamental  mechanistic  modeling  of  the  cellular  calculus  into  a  description  of  population-level  behavior,  and 
thus  to  obtain  a  more  comprehensive  understanding  of  the  immune  system,  with  obvious  implications  for  the 
study  of  disease  detection,  progression,  treatment /control,  etc.  To  that  end,  mathematical  modeling  provides  a 
quantitative  framework  with  which  to  analyze  and  interpret  such  data. 

A  large  number  of  mathematical  models  (see,  e.g.,  the  recent  reviews  [10,  38])  have  been  proposed  with  the 
aim  of  linking  the  generation  structure  (cells  per  number  of  divisions  undergone)  to  quantitative  descriptions  of 
cellular  behavior  (e.g.,  times  to  division  and  death).  Most  recently  [9],  a  class  of  mathematical  models  has  been 
proposed  which  incorporates  the  ‘cyton’  model  [28,  29]  of  cell  division  dynamics  into  a  mathematical  description  of 
flow  cytometry  histogram  data  based  upon  conservation  principles.  Here,  we  revisit  this  new  class  of  models  and 
provide  a  more  complete  discussion  of  some  mathematical  properties  of  the  solutions  which  make  them  amenable 
to  the  fast  computational  approaches  as  described  in  [27].  It  is  also  shown  how  the  new  model  can  be  compared 
with  older  label-structured  models  such  as  those  proposed  in  [12,  27,  42,  47].  Next  the  data  collection  process 
is  considered  in  more  detail  and  a  theoretical  statistical  model  is  derived.  The  mathematical  and  statistical 
models  are  then  incorporated  into  a  rigorous  parameter  estimation  scheme  based  upon  a  weighted  least  squares 
framework  and  members  of  the  proposed  class  of  mathematical  models  are  compared  in  terms  of  their  ability  to 
fit  the  available  data. 


2  CFSE  Data 

CFSE-based  flow  cytometry  experiments  are  performed  by  stimulating  CFSE-labeled  cells  to  divide  by  exposure 
to  either  a  mitogenic  compound  or  a  specific  antigen.  Cells  are  then  placed  into  separate  wells,  one  for  each 
measurement  to  be  made.  Several  protocols  exist  and  can  be  tailored  to  the  needs  of  a  particular  experiment. 
See,  e.g.,  [36,  37,  41,  49,  51].  For  the  experiment  described  here,  peripheral  blood  mononuclear  cells  (PBMCs) 
from  two  healthy  donors  were  stained  with  CFSE  and  stimulated  with  the  mitogen  phytohaemagglutinin  (PHA). 
Measurements  were  carried  out  approximately  every  24  hours  for  five  days,  beginning  one  day  after  stimulation. 

When  a  well  is  selected  for  measurement,  cells  are  additionally  labeled  for  phenotypic  identification  by  anti¬ 
bodies  (anti-CD3  T,  anti-CD4  T  ,  and  anti-CD8  T  cells)  tagged  with  fluorescent  markers.  These  cells  are  then 
analyzed  by  flow  cytometry,  which  records  the  relative  brightness  of  cells  in  various  colors  (corresponding  to 
distinct  fluorescent  markers).  Cells  of  interest  can  then  be  identified  in  the  flow  cytometry  output.  For  this 
experiment,  we  consider  CD4  T  cells  (CD3+,  CD4+,  CD8-)  and  CD8  T  cells  (CD3+,  CD4-,  CD8+).  Once  these 
cells  are  identified,  the  fluorescence  intensity  (in  the  color  channel  consistent  with  CFSE)  is  analyzed  for  each 
cell.  Because  dead  cells  will  disintegrate  shortly  after  death  and  can  be  excluded  by  gating,  mainly  viable  cells 
are  measured  by  the  flow  cytometer  (the  fraction  of  cells  which  are  dead  but  not  yet  disintegrated  is  assumed  to 
be  small). 

The  population  generation  structure  of  the  cells  at  a  given  measurement  time  can  be  visualized  by  organizing 
the  fluorescence  intensity  measurements  of  individual  cells  into  a  histogram,  as  shown  in  Figures  1  and  2.  Because 
of  physical  limitations,  only  a  fraction  of  the  cells  contained  in  a  selected  well  are  actually  analyzed  by  flow 
cytometry.  In  order  to  estimate  the  total  number  of  cells  in  the  measurement  sample,  a  known  number  of 
fluorescent  beads  are  placed  in  each  sample;  these  beads  can  be  identified  and  counted  in  the  flow  cytometry 
output  and  the  ratio  of  beads  counted  to  total  beads  introduced  provides  an  estimate  of  the  fraction  of  the 
sample  acquired  by  the  flow  cytometer.  The  histogram  profiles  obtained  from  the  measured  cells  can  then  be 
normalized  by  the  reciprocal  of  this  quantity. 

It  is  assumed  that  each  well  contains  an  identical  population  of  cells  at  all  times,  that  the  fraction  of  measured 
cells  is  representative  of  the  population  of  cells  in  that  well,  and  that  the  fraction  of  the  total  well  actually  measured 
is  accurately  estimated  by  bead  counting  (though  it  is  possible  to  consider  errors  in  bead  counts;  see  Section  4). 
Under  these  assumptions,  the  total  mass  of  CFSE  should  be  conserved  in  the  histogram  data.  For  this  reason, 
the  mathematical  models  proposed  to  describe  CFSE-based  flow  cytometry  data  (Section  3)  are  derived  from 
standard  conservation  principles  (see,  e.g.,  [13]). 
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Data,  t  =23.5  hrs  Data,  t  =23.5  hrs 


Figure  1:  Histogram  data  for  CD4  T  cells  from  Donor  1  (left)  and  Donor  2  (right).  An  initially  unimodal 
distribution  of  fluorescence  intensity  becomes  multimodal  as  cells  divide  asynchronously.  By  day  5,  subsequent 
generations  of  cells  are  no  longer  detectable  as  fluorescence  resulting  from  CFSE  has  been  diluted  to  the  level  of 
background  autofluorescence. 
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Data,  t  =23.5  hrs 


Data,  t  =23.5  hrs 


Figure  2:  Histogram  data  for  CD8  T  cells  from  Donor  1  (left)  and  Donor  2  (right).  An  initially  unimodal 
distribution  of  fluorescence  intensity  becomes  multimodal  as  cells  divide  asynchronously.  By  day  4,  subsequent 
generations  of  cells  are  no  longer  detectable  as  fluorescence  resulting  from  CFSE  has  been  diluted  to  the  level  of 
background  autofluorescence. 
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A  complete  discussion  of  the  mathematically  relevant  aspects  of  the  data  collection  procedure  can  be  found 
in  [10].  The  goal  of  the  mathematical  modeling  process  is  to  link  a  mathematical  description  of  cellular  division 
and  death  processes  at  the  population  level  to  the  observed  fluorescence  intensity  profiles  as  measured  by  a  flow 
cytometer  (Figures  1  and  2).  Because  each  peak  in  the  flow  cytometry  data  represents  a  cohort  of  cells  having 
completed  the  same  number  of  divisions,  it  is  hypothesized  that  flow  cytometry  data  collected  sequentially  in 
time  for  cells  from  a  single  donor  will  contain  sufficient  information  to  analyze  the  dynamic  response  of  those 
cells  to  stimulus.  This  dynamic  response  can  only  be  accurately  understood  in  the  context  of  a  mathematical 
model  of  the  biological  system,  as  well  as  a  statistical  model  linking  the  mathematical  model  to  the  data.  Such  a 
model  must  be  able  to  account  for  the  slow  natural  loss  of  CFSE  fluorescence  intensity  over  time,  the  dilution  of 
fluorescence  intensity  by  division,  and  the  asynchronous  nature  of  cellular  division  and  death  processes. 


3  Mathematical  Modeling  of  CFSE  Data 

We  begin  by  summarizing  a  partial  differential  equation  model  structured  by  (continuous)  fluorescence  intensity 
and  (discrete)  division  number  which  has  been  proposed  to  describe  histogram  data  from  CFSE-based  proliferation 
assays  [12,  27,  42,  47].  We  then  summarize  a  new  class  of  models  incorporating  cyton  dynamics  into  a  label- 
structured  framework  and  consider  several  different  versions  of  the  cyton  model  at  greater  length.  Finally,  the 
role  of  cellular  autofluorescence  is  briefly  considered. 

3.1  Previous  Label-Structured  Model 

Let  x)  be  the  structured  density  (cells  per  unit  of  fluorescence  intensity)  of  a  cohort  of  cells  having  completed 
i  >  0  divisions  at  time  t  and  with  x  units  of  fluorescence  intensity  resulting  from  induced  CFSE  (that  is,  ignoring 
the  contributions  of  cellular  autofluorescence) .  It  is  assumed  that  this  fluorescence  is  directly  proportional  to  the 
mass  of  CFSE  within  a  cell,  and  thus  can  be  treated  as  a  mass-like  quantity.  These  cells  are  assumed  to  divide 
with  time-dependent  exponential  rate  a.j(f)  and  die  with  time-dependent  exponential  rate  Sif  t).  Then  the  entire 
population  of  cells  can  be  described  by  the  system  of  partial  differential  equations 

=  _  Mt) + 

!  =  _  (Ql(t)  +  m)ni(Ux)  +  Rl{t,x)  (3.1) 


The  recruitment  terms  describe  the  symmetric  division  of  CFSE  upon  mitosis  and  are  given  by  Ri(t,x)  = 
4aj_i(f)?rj_i(f,  2x)  for  i  >  1;  the  form  of  these  recruitment  terms  arises  naturally  from  the  derivation  of  the 
above  system  of  equations  from  conservation  principles  [47] .  The  advection  term  describes  the  rate  of  loss  of  fluo¬ 
rescence  intensity  (resulting  from  the  turnover  of  CFSE),  which  is  assumed  to  depend  linearly  on  the  fluorescence 
intensity  x  with  time-dependent  rate  function  v(t).  This  follows  the  convention  of  [27],  and  includes  exponential 
loss  (v(t)  =  c)  and  Gompertz  decay  (v(t)  =  ce~kt)  as  special  cases.  The  loss  of  CFSE  has  been  observed  to  be 
very  rapid  during  the  first  24  hours  after  initial  labeling  and  much  slower  thereafter  [12,  47]).  Thus  when  data 
is  collected  in  the  first  24  hours,  it  is  more  accurate  [7]  to  describe  the  rate  of  loss  of  fluorescence  intensity  with 
a  time- varying  rate  (e.g.,  Gompertz  decay).  Such  rates  are  consistent  with  the  sequence  of  chemical  reactions 
known  to  occur  during  the  labeling  process  [18].  If  data  is  not  collected  in  the  first  day  after  labeling  with  CFSE 
(as  in  the  data  collected  for  this  report)  then,  as  we  shall  see  below,  exponential  decay  is  sufficient.  For  the 
remainder  of  this  report,  we  assume  v(t)  =  c. 

The  initial  conditions  for  the  model  (3.1)  are 


rii(t0,x) 


$(2:),  i  =  0 
0,  i  >  1  ’ 


(3.2) 


for  x  >  0.  Note  that  a  no-flux  condition  at  x  =  0  is  naturally  satisfied  by  the  form  of  the  advection  term  provided 
n,;(t,  0)  is  finite  (so  that  the  flux  term  is  well-defined)  for  all  i  and  for  all  t  >  0.  The  solution  of  (3.1)  can  be 
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computed  by  the  method  of  characteristics  [47].  Alternatively,  the  following  characterization  of  the  solution  is 
given  in  [42]. 

Proposition  3.1.  The  solution  to  (3.1)  can  be  factored  as 


ni(t,x)  =  Ni(t)rii(t,x). 


The  functions  N(t)  indicate  the  number  of  cells  having  completed  i  divisions  at  time  t  and  satisfy  the  weakly 
coupled  system  of  ordinary  differential  equations 


=  -  (ai(t)  +  ft(i))iVi(i)  +  2ao(t)N0{t) 

=  -  Mf)  +  Pi(t))Ni(t)  +  2ai-i(t)Ni_i(t) 


(3.3) 


with  initial  conditions  No(to)  =  No,  Ni(to)  =  0  for  all  i  >  1.  The  functions  Ui(t,x),  describe  the  distribution  of 
CFSE  within  a  generation  of  cells.  Each  satisfies  the  equation 


dni(t, : 

dt 


-v(t) 


d[xni(t ,  a;)] 
dx 


=  0 


(3.4) 


for  x  >  0  with  initial  condition 


Again,  the  no  flux  condition  at  x 
definition, 


Ui(to,x)  = 


2i$(2ix) 
An  ' 


=  0  is  trivially  satisfied  by  the  form  of  the  advection  term.  Note  that,  by 


No  —  I  <&(x)dx. 
Jo 


We  remark  that  the  form  of  the  equations  presented  above  is  slightly  different  from  that  considered  in  [12] 
as  here  we  initially  neglect  considerations  of  autofluorescence.  The  formulation  above  follows  the  work  of  [27] 
and  allows  for  a  more  intuitive  formulation  of  the  model,  as  well  as  the  fast  numerical  techniques  discussed  below 
and  in  the  appendix.  The  system  above  is  derived  in  terms  of  the  fluorescence  intensity  resulting  only  from 
induced  CFSE ;  the  experimentally  measured  fluorescence  intensity  is  the  sum  of  this  quantity  and  the  cellular 
auto  fluorescence  which  results  from  the  light  absorption  and  emission  properties  of  intracellular  molecules.  Let 
fiift,  x)  be  a  structured  density  for  cells  having  completed  i  divisions  at  time  t,  with  measured  fluorescence  intensity 
x.  While  the  measured  fluorescence  intensity  x  is  given  by  the  sum  of  the  induced  fluorescence  x  and  the  cellular 
autofluorescence,  this  latter  quantity  may  vary  from  cell  to  cell  in  the  population.  As  such,  given  the  solutions 
ni{t,x)  for  *  >  0  to  (3.1),  one  computes  the  densities  hi(t,x)  using  the  convolution  integral  [27,  42] 


hi(t,  x)=  ni(t,  x)p(t,  x  —  x)dx  =  /  x)p(t,  x  —  x)dx , 


(3.5) 


where  p(t ,  £)  is  (for  fixed  time  t)  a  probability  density  function  describing  the  distribution  of  autofluorescence  in 
the  population  (see  [12,  40,  47]  and  Section  3.4  below).  Fast  approximation  techniques  for  the  convolution  (3.5) 
have  been  demonstrated  in  [27]  and  the  references  therein.  These  techniques  are  summarized  in  the  appendix. 
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3.2  New  Class  of  Label-Structured  Models 


While  the  model  (3.1)  for  computing  population  label  structure  has  been  shown  to  accurately  fit  experimental 
data  [12],  it  lacks  a  certain  intuitive  appeal  in  that  the  ‘time-dependent  exponential  rates’  of  division  and  death, 
ctj (t)  and  are  not  explained  in  biologically  relevant  terms  (e.g.,  times  to  division  and  death).  An  alternative 

to  the  mathematical  model  (3.3)  is  the  cyton  model  for  division  dynamics  [28,  29]  which  relates  the  number  of 
cells  in  a  population  directly  to  probability  distributions  describing  times  at  which  cells  divide  or  die.  The  cyton 
model  is  motivated  by  the  assumption  of  independent  regulation  by  the  cellular  machinery  of  times  to  division 
and  death.  Using  the  definition  of  Ni(t)  above,  the  cyton  model  is  described  by  the  set  of  equations 


N0(t)  =  N0 - 


+  nd0'e(s))  ds, 


Ni  (i) 


(3.6) 


where  nflv(t)  and  nfle(t)  indicate  the  rates  (cells  per  hour)  at  which  cells  (having  already  undergone  i  divisions) 
divide  and  die,  respectively,  at  time  t.  For  undivided  cells,  let  0o(i)  and  ipo(t)  be  probability  density  functions 
describing  the  distribution  from  which  the  times  to  division  and  death,  respectively,  are  drawn.  Let  Fq  ,  called  the 
initial  precursor  fraction,  be  the  fraction  of  undivided  cells  which  would  hypothetically  divide  in  the  absence  of 
any  cell  death.  (It  is  assumed  that  in  each  generation,  non-progressing  cells  may  die  according  to  the  probability 
density  function  but  may  not  divide.)  Then,  under  the  assumptions  of  the  cyton  model,  it  follows  that 


nolv(t)  =  F0N0  £  V’o  (s)ds'Sj  (t), 
=  No  (f  -  Fo  Jt  0o (s)ds^  ip0 (*)• 


(3.7) 


Similarly,  one  can  define  probability  density  functions  4>i(t)  and  for  times  to  division  and  death  (measured 
in  hours  since  completion  of  the  ( i  —  1  )th  division/death),  respectively,  for  cells  having  undergone  i  divisions,  as 
well  as  the  progressor  fractions  F)  of  cells  which  would  complete  the  ith  division  in  the  absence  of  cell  death. 
Then  the  cell  division  and  death  rates  are  computed  as 


nfv 


(t)  =  2F.i  f  nf™i(s)  (l  -  f  ipiiOdZ  )  -  s)ds, 

Jt0  \  JO 


ft  /  rt-s 

_ die(j.\  _  o  /  ^div  . 


niie(t)  =  2  /  ni-i(s)  1  -  Fi  /  )  fait  -  s)ds. 

J  to  V  d  0 


(3.8) 


Given  the  success  of  the  cyton  model  in  describing  cell  dynamics,  as  well  as  the  experimental  evidence 
supporting  it  [10,  28,  29],  the  cyton  model  has  been  incorporated  into  a  label-structured  framework  [9]  similar  to 
(3.1).  The  mathematical  ideas  rely  heavily  upon  the  separability  of  the  model  solution  (Proposition  3.1)  originally 
demonstrated  by  [27,  42].  Let  rq(f,  a:)  be  a  structured  density  as  before.  Consider  the  system  of  partial  differential 
equations 


~  V<yt^dx^  =  ~  (noW(*)  +noie(t))  no{t,x) 
dUl  -  =  (2 ndiv(t)  -  nfv{t)  -  nfe(<))  ?h(t,x) 


dt 


(3.9) 


with  initial  conditions  specified  as  for  equation  (3.1).  The  terms  rij(i,  x)  are  described  as  in  (3.4).  This  system  of 
equations  incorporates  the  cyton  model  (3.6)-(3.8)  for  cell  population  dynamics  into  a  label-structured  framework 
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for  use  with  histogram  data.  We  remark  that  (3.9)  corrects  a  typographical  error  in  [9],  where  the  factors  a:) 
were  missing  from  the  right  side  of  the  equations. 

Proposition  3.2.  The  solution  of  (3.9)  is 

rii(t,  x)  =  Ni{t)hi{t,  x) 

where  the  quantities  Ni(t)  satisfy  (3.6)  and  fii(t,x)  satisfy  (3.4). 

Proof.  The  proof  follows  immediately  by  the  direct  substitution  of  the  stated  solution  into  (3.9).  Working  with 
the  left  side  of  (3.9)  for  the  ith  equation, 

dnj{t,x)  u^d[xn.i(t,x)\  _  d[Ni(t)fii(t,x)]  _t^d[xNi(t)fii(t,x)] 
dt  dx  dt  dx 

dNi{t)  f  dni(t,  x)  d[xni(t,x)\ 

=  -nrn'{t’x)  +  Ndt)  (— s—  -  ”(t)  a, 

=  (2n?i”1-nf”-nf')sJ(t,x), 

which  is  exactly  the  right  side  of  (3.9).  For  the  purposes  of  this  proof,  it  is  assumed  that  ri:lf'{  =  0  so  that  the 
equations  for  no(t,  x)  are  well-defined.  It  is  easy  to  check  that  the  initial  and  boundary  conditions  for  (3.9)  are 
satisfied  by  the  above  solution.  □ 


Given  the  densities  rii(t,x)  computed  according  to  (3.9),  one  can  compute  the  densities  hi(t,x)  via  the 
convolution  (3.5)  as  before.  Much  like  the  original  model  (3.1),  the  new  model  (3.9)  can  be  fit  directly  to 
histogram  data  from  CFSE-based  experiments  and  is  highly  accurate  [9] .  Significantly,  the  new  model  describes 
the  dynamics  of  a  dividing  population  of  cells  in  intuitive  terms  (i.e.,  probability  distributions  of  times  to  divide 
and  die).  This  is  the  primary  advantage  of  the  new  class  of  models  over  the  previous  modeling  framework. 
Similarly,  while  the  cyton  model  has  been  widely  used  to  analyze  cell  count  data  obtained  from  CFSE  data  (e.g., 
through  a  deconvolution  process;  see  [10]),  the  new  class  of  models  can  be  fit  directly  to  CFSE  histogram  data. 
As  a  result,  the  class  of  models  is  less  dependent  upon  peak  separation  or  a  high  frequency  of  cells  which  respond 
to  stimulus.  Moreover,  the  fit  of  the  model  to  data  can  be  assessed  in  a  statistically  rigorous  manner  (see  Section 
4  below). 

Although  the  motivation  for  this  model  formulation  is  clear  (combining  cyton  and  label  dynamics  in  a  division- 
dependent  compartmental  model)  the  form  of  the  new  model  is  complex,  describing  the  population  densities 
n,;(t,  a:)  in  terms  of  yet  another  set  of  density  functions,  iii(t,  x).  A  simple  reformulation  shows  that  the  new  class 
of  models  is  consistent  with  the  mass-conservation  principles  of  the  old  label-structured  model.  Moreover,  this 
reformulation  shows  how  the  two  model  forms  can  be  related  and  directly  compared. 

Recall  the  definitions  of  n.i(t,x),  Ni(t),  and  hi(t,x)  given  in  Proposition  3.2.  Note  that 


x)dx  =  Ni(t), 


and 


n(t,  x)dx  =  / 

Jo 


Ui(t,x) 


dx  = 


1 


rii(t,  x)dx  =  1. 


io  Jo  Ni(t)  Ni{t)  J o 

Thus  the  quantities  hi(t,x)  can  be  considered  as  probability  density  functions  for  the  distribution  of  CFSE  in 
cells  having  divided  i  times  at  time  t.  When  considering  cell  death,  the  mathematical  terms  nfle(t)fii(t,  x)  in 
(3.9)  reflect  the  tacit  assumption  that  the  rate  at  which  cells  die  ( nfie(t ))  is  independent  of  the  label  distribution 
hi(t,x)  of  those  cells.  Moreover, 


nfe(t)m(t,x)  =  nfe(t) 


Ui(t,x) 

Ni(t) 

-  jv,(t) 

=  /3i(t)rii(t,x), 


-.die 


W 


with 
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Pi{t) 


»fe(0 

Ni(t)  ’ 


which  is  exactly  the  same  form  as  (3.1).  Similar  statements  hold  true  for  the  rates  of  dividing  cells  and  the  terms 
ctj (t)  in  Equation  (3.1)  if  Ft  of  (3.7)-(3.8)  equal  1  for  all  i.  For  0  <  Fi  <  1,  then  this  is  a  more  complex  issue 
and  indeed  is  the  subject  of  some  of  our  current  efforts.  This  fact  can  be  used  as  the  basis  for  a  quantitative 
comparison  of  the  two  model  formulations,  as  well  as  to  give  physical/biological  meaning  to  the  time  dependent 
exponential  rates  a.i(t)  and  Pi(t)  used  previously  to  describe  cell  division  and  death. 


3.3  The  Cyton  Class  of  Models 

It  follows  from  the  form  of  the  equations  (3.6)-(3.9)  that  the  generation  structure  of  the  population  (cells  per 
division  number)  is  completely  determined  by  the  functions  <j>i(t)  and  'ipi(t)  and  the  progressor  fractions  f7).  To 
motivate  the  form  of  the  functions  4>i{t)  and  define  the  random  variable  Tdm  to  be  the  time  required  for  a 

progressing  cell  to  complete  the  ith  division,  with  the  clock  starting  from  the  completion  of  the  (i  —  l)th  division. 
(That  is,  the  random  variables  Tdm  are  defined  in  the  temporal  reference  frames  of  the  individual  cells.)  Similarly, 
define  the  random  variables  Tfle  to  be  the  time  required  for  a  newly  divided  cell  to  die.  The  cyton  model  is 
built  from  the  premise  that  these  two  random  variables  are  independent.  Upon  the  completion  of  the  ith  division 
(or  upon  activation,  for  i  =  0),  every  cell  realizes  a  new  value  for  Tdlv  and  T*e;  whichever  realization  is  smaller 
determines  the  eventual  fate  of  the  cell.  The  functions  </>,:(£)  and  ipiif)  are  the  probability  density  functions  for 
Tflv  and  Tdte ,  respectively  (which  are  assumed  to  be  common  for  all  cells  having  completed  i  divisions). 

Experimental  evidence  suggests  that  the  functions  and  'fi(t)  can  be  heuristically  described  by  lognormal 
probability  density  functions  [28,  29].  Thus  for  all  t  >  0, 


todivV 2^CXP  ( 
tafeV 2)feXP  ( 


(log  t  Pi™ ) 2  \ 

2 (erf™)2  )  ’ 

(log  t  -  M?e)2\ 

2(^fe)2  )  ’ 


(3.11) 


where  the  parameters  /zf™  and  erf™  represent  the  means  and  standard  deviations  of  the  natural  logarithms  of 
the  random  variables  Tdlv  (and  similarly  for  Tdle).  Since  it  is  more  intuitive  to  discuss  the  means  and  standard 
deviations  of  the  random  variables  Tfm  and  Tdle  directly  (as  opposed  to  the  means  and  standard  deviations  of 
their  logarithms)  these  quantities  are  easily  defined  in  terms  of  the  parameters  /zf™,  pdle,  erf™,  and  <rf*e: 


/  (ndiv)2\ 

E[T? iv]  =  exp  [p?v  + 

/  (rrdie)2\ 

E[Tdie]  =  exp  (m?‘ 3  + 

Var[T?iv }  =  (exp  ((af™)2)  -  l)  exp  (2pf™  +  (af™)2) 
Var[Tdie]  =  (exp  ((afe)2)  -  l)  exp  (2/zf e  +  (af e)2)  . 


For  the  basic  cyton  model,  it  is  standard  (following  the  work  [28,  29])  to  assume  that  the  random  variables 
T“™  are  identically  distributed  for  all  i  >  1  and  that  the  random  variables  Tdle  are  identically  distributed  for  all 
i  >  1.  These  distributions  may  be  different  from  the  corresponding  random  variables  for  undivided  cells  (i  =  0). 
Thus 
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It  is  also  assumed  that  fq  =  1  for  all  i  >  1  in  the  basic  cyton  model. 

Of  course,  any  number  of  generalizations  of  the  basic  cyton  model  is  possible.  For  instance,  following  [28],  the 
fractions  Fi  can  be  defined  in  terms  of  a  division  destiny.  Among  the  cells  which  are  activated  to  divide  (FqNq 
of  them),  let  pi  be  the  probability  that  a  cell  (or  its  progeny)  ceases  to  be  activated  after  completing  i  divisions 
and  define  the  cumulative  probabilities 

i 

=  ^PJ- 

3=  1 

(Note  that  we  must  have  Cj  — »•  1  as  i  — >  oo.)  It  follows  that  the  progressor  fractions  (for  i  >  1)  are 


Fi  = 


l  -a 

1—Ci-l  ’ 


0, 


Ci-1  <  1 
Ci-1  =  1. 


(3.12) 


Rather  than  estimate  the  progressor  fractions  Fj  (or  the  probabilities  pf)  independently,  we  follow  the  approach 
suggested  in  [28]  and  assume  that  the  probabilities  pi  can  be  described  as  a  discrete  normal  density  function 
defined  on  the  nonnegative  integers.  Thus  the  values  of  the  probabilities  pi  (and  the  progressor  fractions  Fi) 
are  uniquely  determined  by  the  mean  D ^  and  the  standard  deviation  Da  of  a  discrete  normal  distribution.  This 
assumption  has  been  shown  to  be  consistent  with  experimental  data  [28]  and  has  the  beneficial  effect  of  reducing 
the  total  number  of  parameters  of  the  mathematical  model. 

We  can  now  define  the  division  destinies  in  terms  of  the  progressor  fractions  (3.12).  The  division  destiny 
di  is  defined  to  be  the  fraction  of  cells  out  of  those  cells  in  the  original  population  which  would  have  proceeded 
through  exactly  i  divisions  in  the  absence  of  any  cell  death.  These  quantities  are  computed  as 


di  = 


1  ~  F0,  i  =  0 
F0pi,  i  >  1. 


(3.13) 


It  should  be  noted  that  this  definition  does  not  make  any  assumptions  regarding  the  exact  lineage  of  cells  (which 
cannot  be  determined  from  flow  cytometry  data)  so  that  the  progeny  of  a  single  cell  are  not  assumed  to  all  undergo 
the  same  number  of  divisions.  Rather,  one  can  consider  a  fractional  number  of  precursors.  For  example,  consider 
a  single  cell  which  divides  once,  after  which  one  of  the  two  daughter  cells  divides  again.  Then  there  are  three 
cells  in  the  total  population,  and  the  division  destinies  are  do  =  0,  d\  =  1/3,  d2  =  2/3.  Though  counterintuitive 
for  a  single  cell,  division  destinies  provide  an  indication  of  the  number  of  divisions  undergone  averaged  over  the 
population  of  precursors. 

Following  the  work  presented  in  [9],  one  may  also  generalize  the  death  rate  mechanism  for  undivided  cells 
to  incorporate  a  separate  set  of  behaviors  for  unactivated  cells.  In  particular,  it  can  be  assumed  that  a  fraction 
Pd  of  such  cells  will  remain  dormant  and  neither  divide  nor  die  during  the  experiment.  The  remaining  fraction 
(1  —  pd)  will  die  with  some  exponential  rate  /?  which  is  independent  of  the  death-rate  distribution  of  activated 
cells.  It  follows  that  the  probability  density  function  describing  cell  death  for  undivided  cells  is 


i’o  (t) 


Fn 


ta^V^K 


exp 


(log  l  -  H oie)2\ 

2  «e)2  ) 


+  (l-pd)(l~Fo)pe~0t. 


(3.14) 


It  should  be  noted  that  this  generalization  changes  the  interpretation  of  the  random  variable  T*e  and  its  rela¬ 
tionship  to  the  parameters  and  &ole  in  the  sense  that  the  parameters  pole  and  cr*6  describe  the  statistical 
properties  of  progressing  cells  only. 

While  the  more  complex  death  rate  function  (3.14)  was  found  to  accurately  describe  a  CFSE  data  set  in  [9], 
the  data  sets  collected  for  this  manuscript  differ  in  that  the  first  measurement  was  taken  approximately  24  hours 
after  stimulation  by  PHA  (as  opposed  to  immediately  following  stimulation).  As  a  result,  the  initial  condition  for 
the  mathematical  model  represents  only  those  cells  which  have  not  died  in  the  first  24  hours  after  stimulation. 
It  seems  reasonable  to  hypothesize  that  such  cells  are  unlikely  to  die  at  subsequent  measurement  times.  Thus 
an  additional  possibility  is  ifoit)  =  0  for  all  t.  (This  ipo(t)  is  not  a  proper  probability  density  function,  but  is 
sufficient  to  describe  the  intended  behavior.  Equivalently,  one  could  assume  the  function  'fo(t)  has  a  large  mean 
and  small  variance  so  that  the  support  of  the  density  function  is  effectively  limited  to  a  region  beyond  the  final 
measurement  time.) 
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Model 

Description 

Parameters  (cyton  only) 

Model  1 

Basic  cyton  model;  Equations  (3.11),  Fi  —  1  for  all  *  >  1 

9 

Model  2 

Basic  cyton  model  plus  division  destiny  according  to  (3.12) 

11 

Model  3 

Basic  cyton  model,  but  with  equation  (3.14)  for  undivided  cell  death 

11 

Model  4 

Basic  cyton  model,  but  with  no  undivided  cell  death  =  0) 

7 

Model  5 

Combine  models  2  and  3 

13 

Model  6 

Combine  models  2  and  4 

9 

Model  7 

Model  1,  but  with  equation  (3.15)  for  undivided  cell  division 

12 

Model  8 

Model  2,  but  with  equation  (3.15)  for  undivided  cell  division 

14 

Model  9 

Model  3,  but  with  equation  (3.15)  for  undivided  cell  division 

14 

Model  10 

Model  4,  but  with  equation  (3.15)  for  undivided  cell  division 

10 

Model  11 

Model  5,  but  with  equation  (3.15)  for  undivided  cell  division 

16 

Model  12 

Model  6,  but  with  equation  (3.15)  for  undivided  cell  division 

12 

Table  1:  List  of  the  possible  cyton  model  parameterizations  considered  in  this  report.  These  models  are  compared 
(in  terms  of  their  ability  to  describe  experimental  data  sets)  in  Tables  3  and  4. 


Finally,  we  consider  one  possible  generalization  of  the  density  function  for  time-to- first-division  for  progressing 
cells.  If  there  are  multiple  subpopulations  (e.g.,  naive  vs.  memory  cells)  contained  within  the  population  under 
study,  the  density  <f>o{t)  may  be  multimodal.  For  simplicity,  we  consider  a  bimodal  density  function  which  is  a 
weighted  sum  of  two  lognormal  distributions, 


f  (  (log 

t<v^eXP  V  2«)2  J 


1  -/  (  (log  t-^)2\ 

K^exp  V  2(<b)2  )  ’ 


where  /  £  [0,1]  is  a  weighting  parameter. 

The  possible  model  parameterizations  considered  in  this  report  are  summarized  in  Table  1. 


(3.15) 


3.4  Distribution  of  Cellular  Autofluorescence 

To  this  point  we  have  considered  a  class  of  models  based  on  the  cyton  modeling  framework  which  describe  the 
dynamic  population  generation  structure  for  dividing  cells.  This  class  of  models  has  been  incorporated  into 
a  label-structured  partial  differential  equation  model  derived  by  considering  the  CFSE  in  a  mass-conservation 
framework.  As  discussed  previously,  once  the  structured  densities  riiit ,  x)  (in  terms  of  the  fluorescence  x  resulting 
from  CFSE)  have  been  constructed,  these  quantities  must  be  related  to  the  measured  fluorescence  intensity  x 
(which  includes  the  contribution  of  cellular  autofluorescence)  via  the  convolution  (3.5).  Autofluorescence  is  the 
result  of  the  absorption  and  emission  properties  of  molecules  which  are  naturally  found  within  all  cells  and  is 
present  even  in  the  absence  of  an  added  fluorescent  label.  Mean  autofluorescence  is  known  to  increase  as  cells  are 
activated  to  divide  [1],  probably  as  a  result  of  the  production  of  additional  intracellular  components  associated  with 
increased  metabolic  activity  within  the  cell.  Thus  the  notation  of  (3.5)  explicitly  includes  the  time-dependence 
of  the  autofluorescence  density  function,  p(t,£).  Because  these  intracellular  molecules  are  partitioned  among 
daughter  cells  during  cell  division,  the  distribution  of  autofluorescence  can  be  intuitively  considered  as  a  growth 
and  fragmentation  process,  which  is  known  to  produce  skew-right  density  functions  such  as  the  lognormal  density 
function  [26].  In  fact,  it  has  been  shown  [12]  that  the  distribution  of  autofluorescence  in  the  population  can  be 
well-approximated  using  a  lognormal  density  function,  and  thus  can  be  characterized  by  its  mean  and  its  variance. 
This  observation  has  been  used  as  the  basis  for  approximation  techniques  to  the  convolution  (3.5)  [27]. 

To  test  the  assumption  of  lognormality,  a  portion  of  PBMC  cells  from  each  donor  were  set  aside  and  stimulated 
with  PHA  but  never  labeled  with  CFSE.  Thus  the  measured  fluorescence  distribution  of  these  cells  (represented 
in  histogram  form)  can  be  used  to  approximate  the  density  function  p(t,£)  representing  the  actual  population 
distribution  of  autofluorescence.  This  autofluorescence  data  is  depicted  for  the  two  donors  and  cell  types  in 
Figures  3  and  4.  To  assess  the  lognormal  approximation,  we  used  two  parameter  estimation  schemes  to  construct 
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lognormal  density  functions  from  the  autofluorescence  data.  The  first  method  is  the  method  of  moments,  in  which 
the  exact  mean  and  variance  of  the  measured  cells  was  computed  and  a  lognormal  curve  was  constructed  with 
the  same  mean  and  variance.  (In  the  figures,  the  resulting  lognormal  density  function  has  been  scaled  by  the 
number  of  cells  in  the  data  to  facilitate  comparison.)  The  second  method  is  to  use  least  squares  to  estimate  the 
scale,  mean,  and  variance  of  a  lognormal  density  function. 

Though  the  autofluorescence  data  is  not  perfectly  lognormal,  the  assumption  is  fairly  accurate  for  both  cell 
types  and  both  donors.  For  CD4  T  cells,  the  lognormal  approximation  becomes  more  accurate  as  time  progresses. 
This  is  consistent  with  the  existence  of  an  initial  transient  distribution  of  autofluorescence  which  corresponds  to 
the  quiescent  state  of  the  cells  at  the  start  of  the  experiment;  as  cells  are  activated  to  divide,  the  lognormal  (or 
at  least  skew-right)  distribution  emerges  possibly  as  a  result  of  growth  and  division  processes  [26].  For  CD8  T 
cells,  the  validity  of  the  approximation  does  not  change  much  in  time. 

Significantly,  the  statistical  moments  of  the  autofluorescence  distribution  are  observed  to  change  in  time 
(Figure  5).  This  is  particularly  true  for  the  mean.  The  significant  increase  in  mean  autofluorescence  between 
24  and  48  hours  is  known  to  be  associated  with  cellular  activation.  Though  the  increase  is  large,  it  is  of  little 
consequence  for  mathematical  modeling  because  the  contribution  of  autofluorescence  to  the  fluorescence  intensity 
measurements  is  very  small  (less  than  1%)  for  undivided  cells.  The  decrease  in  autofluorescence  as  measured 
at  approximately  t  =  96  hours  is  more  problematic  as  some  cells  will  have  completed  multiple  divisions  by  that 
time;  the  cause  of  the  decrease  is  unknown.  It  is  possible  that  a  change  in  instrument  settings  between  the  two 
measurement  days  could  account  for  this  effect.  Yet  this  seems  unlikely  as  comparable  changes  are  not  observed 
in  the  data  collected  for  cells  labeled  with  CFSE  (Figures  1  and  2).  Another  possibility  is  nutrient  depletion;  by 
the  third  day  of  the  experiment,  the  cells  begin  to  run  out  of  the  nutrients  originally  placed  in  culture  and  these 
must  be  replaced.  Nutrient  depletion  could  cause  activated  cells  to  die  or  return  to  a  quiescent  state. 

Based  upon  these  observations,  the  population  distribution  of  autofluorescence  can  be  accurately  described 
as  a  lognormal  density  function  at  each  measurement  time.  Thus 


P{t,Q 


1  cxp  (  Ogg 

£rIo(i)v^F  V  2  aXa(t?  )■ 


(3.16) 


where  the  parameters  pXa(t)  and  crXa(t)  are  related  to  the  mean  and  standard  deviation  of  the  autofluorescence 
distribution  (at  time  t)  by  the  formulas 


fe0(t)  =log(i?[a;0(£)])  -  -  log  (  1  + 


^x0W  =log  ^1  + 


STD{xa(t))‘ 

E[xa(t)}2 


STD{xa(t))2  \ 
E[xa(t)}2  ) 


In  the  context  of  parameter  estimation  for  the  mathematical  model  (3.9),  we  consider  two  frameworks. 
The  first  is  to  use  the  method  of  moments  (as  above)  with  the  autofluorescence  data  to  compute  E[xa(t)\  and 
ST D(xa(t))',  these  values  will  then  be  considered  fixed  while  the  remaining  parameters  of  the  mathematical  model 
are  determined  by  using  least  squares  to  fit  the  data  in  Figures  1  and  2  (see  Section  4  below).  The  second  method 
is  to  ignore  the  time-dependence  of  the  distribution  and  assume  E[xa(t)\  =  E[xa],  STD(xa(t ))  =  STD{xa).  The 
values  of  E[xa]  and  STD(xa )  are  estimated  in  a  least  squares  framework  as  described  in  Section  4,  and  these 
estimated  values  can  be  compared  to  the  values  returned  by  the  method  of  moments.  One  could  also  attempt  to 
parameterize  and  estimate  a  function  p(t ,  £)  which  is  time-dependent.  This  is  a  more  complex  estimation  problem 
and  is  unlikely  to  be  identifiable;  therefore  we  do  not  consider  it  here. 
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Figure  3:  Measured  autofluorescence  distributions  in  comparison  to  fitted  lognormal  curves  for  Donor  1  (left)  and 
Donor  2  (right)  for  CD4  T  cells.  The  least  squares  fit  (LS)  to  data  is  visually  better  than  the  fit  obtained  by 
the  method  of  moments  (MM),  though  the  difference  is  small  and  becomes  less  noticeable  at  later  measurement 
times.  This  pattern  is  more  noticeable  for  CD4  T  cells  (here)  than  for  CD8  T  cells  (Figure  4). 
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Figure  4:  Measured  autofluorescence  distributions  in  comparison  to  fitted  lognormal  curves  for  Donor  1  (left)  and 
Donor  2  (right)  for  CD8  T  cells.  The  least  squares  fit  (LS)  to  data  is  visually  better  than  the  fit  obtained  by  the 
method  of  moments  (MM),  though  the  difference  is  small. 
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AutoFI  Mean  CD4  Vivid  AutoFI  Mean  CD8  Vivid 


Figure  5:  Top:  Mean  autofluorescence  as  estimated  by  the  method  of  moments  (solid  lines)  and  least  squares 
(dashed  lines).  Bottom:  Standard  deviation  as  estimated  by  the  method  of  moments  (solid  lines)  and  least  squares 
(dashed  lines).  Data  shown  for  CD4  T  cells  (left)  and  CD8  T  cells  (right).  Mean  autofluorescence  is  observed  to 
increase  significantly  in  the  first  two  days  of  the  experiment,  an  effect  which  is  known  to  be  the  result  of  cellular 
activation.  The  cause  of  the  decrease  in  mean  autofluorescence  at  approximately  96  hours  is  unknown;  nutrient 
depletion  is  one  hypothesized  cause. 
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4  Statistical  Modeling  of  CFSE  Data 

Models  similar  to  those  discussed  in  Section  3  have  previously  been  shown  to  accurately  describe  population 
generation  structure  as  measured  by  flow  cytometry  [9].  However  the  statistical  properties  of  the  measurement 
process  have  not  been  nearly  as  carefully  considered.  A  major  limitation  of  the  current  modeling  framework  lies 
not  in  the  mathematical  model  itself  but  rather  in  the  statistical  model  which  links  the  mathematical  model  to 
the  data.  An  accurate  statistical  model  is  of  vital  importance  for  the  consistent  estimation  of  model  parameters, 
as  well  as  the  unbiased  estimation  of  confidence  intervals  around  those  parameters  [2,  13,  19,  45].  Additionally, 
an  accurate  statistical  model  is  necessary  for  the  rigorous  comparison  of  different  model  parameterizations  and 
generalizations  [3,  5,  16]  and  the  optimal  design  of  experiments  [4]. 

To  this  point,  we  have  discussed  a  class  of  mathematical  models  which  combine  the  cyton  modeling  framework 
of  [28,  29]  with  the  label  and  division  structured  population  models  of  [12,  27,  42,  47]  to  describe  CFSE  data. 
Given  several  members  of  this  class  of  models  (see,  e.g.,  Table  1)  there  is  a  need  to  compare  the  mathematical 
models  on  a  quantitative  basis  in  order  to  identify  which  model  provides  the  ‘best’  (in  an  appropriate  sense) 
description  of  an  underlying  experimental  data  set.  Several  techniques  based  upon  information  theory  [16]  or 
asymptotic  properties  of  least  squares  estimators  [3,  5,  24]  have  been  developed  for  this  purpose.  In  all  cases,  the 
techniques  are  premised  upon  an  accurate  statistical  model  which  links  the  mathematical  model  to  the  collected 
data. 

Mathematical  descriptions  of  CFSE  data  have  generally  described  numbers  of  cells  per  generation  as  estimated 
from  histogram  data,  rather  than  the  histogram  data  itself.  As  such,  little  consideration  has  been  given  to  the 
statistical  model  which  generates  the  histogram  data.  In  their  likelihood  estimation  framework,  Hyrien  and  Zand 
[31]  propose  that  the  marginal  probability  density  of  each  datum  is  normally  distributed,  and  that  the  variance 
of  this  normal  distribution  is  constant  for  all  data  points.  Least  squares  estimators,  though  not  restricted  to  any 
parametric  class  of  probability  density  functions,  have  also  generally  assumed  a  constant  variance  error  model 
[7,  8,  9,  12,  34,  35,  47].  However,  it  has  been  shown  that  such  an  assumption  does  not  accurately  describe  the 
variance  as  observed  in  actual  data  sets  [8,  12,  47].  Another  common  error  model  for  least  squares  estimation, 
in  which  the  variance  of  each  data  point  is  assumed  to  be  directly  proportional  to  the  square  of  the  model  value 
at  that  point  (a  constant  coefficient  of  variance  model)  has  also  been  hypothesized,  but  was  again  observed  to 
be  inaccurate  [8,  12,  47].  Here,  we  revisit  the  discussion  of  [47,  Ch.  4]  to  consider  the  probabilistic  aspects  of 
the  actual  experimental  process  itself  and  derive  a  hypothetical  statistical  model  from  a  theoretical  basis.  This 
statistical  model  is  then  incorporated  into  a  weighted  least  squares  estimation  scheme  and  several  computational 
algorithms  are  proposed. 

4.1  Theoretical  Statistical  Model 

Define  the  structured  densities  hi(t,x)  (in  terms  of  measured  fluorescence  intensity)  for  cells  having  completed  i 
divisions  as  in  Section  3.  Then  the  structured  density  for  the  entire  population  of  cells  is 

n(t,  x)  =  hi(t,  x).  (4.1) 

i 

Because  CFSE  histogram  data  are  most  commonly  represented  using  a  base  10  logarithmic  scale,  we  define  the 
change  of  variables  z  =  log10(i)  to  arrive  at 

n(t,  z)  =  10z  log(10)n(t,  10z) 

which  gives  the  structured  density  (measured  in  cells  per  base  10  log  unit  intensity)  for  the  entire  population 
of  cells  at  time  t.  Let  Nk  be  a  random  variable  representing  the  number  of  cells  measured  at  time  tj  with  log- 
fluorescence  intensity  in  the  region  [zk,Zk+ 1).  The  goal  of  the  statistical  model  is  to  link  the  structured  density 
h(t,  z)  to  the  data  random  variables  Nl  in  a  manner  that  is  consistent  with  the  statistical  properties  of  the  random 
variables  Nl.  Moreover,  in  the  experimental  process,  the  quantities  Nl  are  not  directly  measured.  Rather,  only 
a  fraction  of  the  contents  of  each  measurement  well  are  acquired  by  the  measurement  apparatus.  In  order  to 
account  for  this  discrepancy,  a  known  number  of  beads  (which  can  be  identified  and  counted  in  the  flow  cytometry 
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output)  are  contained  within  each  sample  to  be  measured.  By  comparing  the  number  of  beads  acquired  to  the 
number  of  beads  known  to  be  in  the  tube,  one  is  able  to  estimate  the  fraction  of  the  sample  acquired. 

Let  q  be  the  vector  of  parameters  of  the  mathematical  model  (that  is,  q  contains  the  parameters  necessary 
to  describe  the  cyton  dynamics  as  well  as  the  parameters  describing  the  label  loss  function  v(t)  and  the  autofluo¬ 
rescence  distribution  p(t,  £))  so  that  we  may  rewrite  n(t ,  z)  =  h(t,  z;  q).  In  order  to  derive  an  error  model  for  the 
histogram  data  we  first  make  the  common  assumption  that  the  model  is  correctly  specified  so  that  the  structured 
population  density  h(t,z;qo)  (where  <fo  is  a  hypothetical  ‘true’  parameter)  perfectly  describes  the  population  of 
cells.  Let  iVj(f)  be  defined  as  in  (3.6)  and  let  N(t)  =  ^iVi(t).  That  is,  N(t)  is  the  total  number  of  cells  in  the 
population  at  time  t.  Define 

_  n{tJlZ-q0) 

M>~  JV(tj)  ■ 

It  follows  that  Pj(z)  is  a  probability  density  function.  Let  Sj  be  the  number  of  cells  of  interest  (e.g.,  CD4  T  cells) 
sampled  at  measurement  time  tj.  Then  one  can  consider  the  sample  of  Sj  cells  (of  interest)  to  be  taken  without 
replacement  from  the  total  population  of  N(tj)  cells;  the  fluorescence  intensity  of  the  sampled  cells  is  subject  to 
the  sampling  density  pj(z).  It  should  be  carefully  noted  that  there  are  numerous  steps  required  to  separate  the 
cells  of  interest  from  the  actual  culture  of  cells  passing  through  the  cytometer.  See,  e.g.,  [10,  Sec.  2],  References 
to  the  total  number  of  cells  N(t),  and  the  number  of  sampled  cells  Sj  are  understood  to  refer  only  to  the  specific 
cells  of  interest  in  the  experiment.  For  the  moment,  we  make  the  additional  assumption  that  these  two  numbers 
are  exact  and  are  not  subject  to  any  errors  (systematic,  experimental,  or  otherwise)  caused  by  gating,  etc. 

Let  B  be  the  total  number  of  beads  (in  each  sample  tube)  which  are  used  to  quantify  the  fraction  of  the 
population  of  cells  which  is  measured  at  time  tj  and  let  bj  be  the  ‘true’  number  of  beads  passing  through  the 
cytometer.  By  this,  we  mean  the  exact  number  of  beads  which  would  pass  through  the  cytometer  if  the  measured 
culture  were  perfectly  homogeneous,  etc.  It  follows  that 


Sj  = 


Now  consider  the  kth  histogram  bin  [zk,zk+ 1).  The  number  of  cells  in  the  whole  population  which  are 
contained  in  this  bin  is 

rZk+ 1 

(4.2) 


rzk+ 1 

I[h](tj,zk-,qo)  =  /  h(tj,z;q0)dz. 

J  Zh. 


Let  Mk  be  a  random  variable  representing  the  number  of  cells  (out  of  the  sampled  population)  counted  into 
the  kth  bin.  Because  the  measurement  process  represents  a  sampling  without  replacement,  it  follows  that  Mk  is 
described  by  a  hypergeometric  distribution, 


Mk  ~  HypG(N(tj),I[n\(tj,zk;q0),Sj). 


That  is,  Sj  cells  are  sampled  without  replacement  from  a  population  containing  a  total  of  N(tj)  cells,  of  which 
I[n](tj,Zk',qo)  are  of  interest.  We  make  the  following  assumptions  regarding  the  measurement  process: 

•  N(tj)  »  Sj  (and  thus  I[n(tj,  zk;qo)]  » 


U  <  6  S  N (tj) 


<  1  -  e  <  1. 


Then  it  can  be  shown  [23]  that  MJk  dlstbny  M3k  where 


A~>j  Kr  ( Sjl[n]{tj,zk\q0)  Sjl[h]{tj,zk\q0)  f  _  I[n](tj,zk-,q0)\  \ 
fc  V  N(tj)  ’  N(tj)  iv  N(tj)  )) 

-  \r  fbj  N(tj)Iifl\(tj^zk;qo)  bj_  N(tj)I[n}(tj,zk;  q0)  /  I[n]{tj,zk;q0)\  \ 

\B  N(tj)  B  N{tj)  V  N(tj)  JJ 
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Notation 

Description 

h(t,z) 

Log-transformed  label  structured  density 

N(t) 

Total  number  of  cells  in  the  population  at  time  t 

pAz) 

Probability  density  function  from  which  cells  are  sampled 

Sj 

Number  of  cells  sampled  at  time  tj 

B 

Total  number  of  beads  originally  placed  into  each  well 

bi 

‘True’  number  of  beads  counted  by  the  cytometer  at  time  tj 

f  O 

1 — 1 

‘True’  number  of  cells  from  the  total  population  belonging  in  the  kth  histogram  bin  at  time 

Aj 

Random  variable  representing  the  number  of  cells  counted  into  the  kth  histogram  bin  at 
time  tj 

K 

bj 

Actual  number  of  beads  counted  by  the  flow  cytometer  (a  realization  of  bj) 

K 

Random  variable  resulting  when  Ml  is  scaled  by  the  ratio  B /bj 

K 

The  actual  data,  a  realization  of  IV- 

A, 

Random  variable  representing  the  bead  count  error  ratio  bj/bj 

Table  2:  Summary  of  notation  for  the  statistical  model. 


The  final  approximation  is  valid  provided  I[h](tj,  zk;qo)  «  N(tj),  which  is  a  perfectly  reasonable  assumption. 

It  can  easily  be  shown  that  the  first  assumption  regarding  the  measurement  process  is  accurate  provided 
bj/B  <<  1,  which  again  is  reasonable  (the  ratio  is  typically  less  than  0.1).  This  assumption  is  necessary  to 
ensure  that  the  sampling  (without  replacement)  process  is  conducted  in  a  such  a  way  that  the  ratio  of  cells  of 
interest  to  total  cells  is  approximately  constant  during  the  measurement.  The  assumption  is  unusual  in  that 
it  places  a  restriction  on  the  total  amount  of  data  which  can  be  collected.  The  second  assumption  regarding 
the  measurement  process  bounds  the  probability  that  a  cell  belongs  to  a  particular  bin  away  from  zero  and  one 
(although  this  assumption  is  not  strictly  necessary  in  some  cases  [33]).  In  practice,  this  assumption  is  only  violated 
when  I[h(tj,zk;  <fo)]  «  0. 

Finally,  when  the  measurements  are  actually  taken,  a  certain  number  of  beads  bj  are  actually  counted.  We 
would  certainly  hope  that  bj  ss  bj ;  however,  we  can  think  of  bj  as  a  realization  of  some  random  variable  (which  may 
or  may  not  be  an  unbiased  estimator  of  bj,  depending  upon  any  systematic  error  that  might  occur  in  obtaining 
bead  counts  from  flow  cytometry  data).  To  obtain  the  histogram  data  which  is  actually  used  to  calibrate  the 
mathematical  model,  one  scales  the  sampled  cell  counts  Ml  by  the  inverse  of  the  fraction  of  the  total  population 

number  of  counted  beads  bj).  Thus  the  data  may  be  represented  by  the 

K  =  ymI 

bj 

(B  b  B ^  b  ■  \ 

y  ^HnKtj^k;  qo),  -^^I[n}(tj,zk;  q0)j 

I  Aj/[h]  (tj ,  zkj  qo),  Aj j— 

V  bJ 

where  we  have  defined  A j  =  bj/bj.  The  quantities  A j  in  effect  represent  a  ‘scaling  error’  and  are  sufficient  to 
explain  the  apparent  violation  of  conservation  principles  often  noticed  in  flow  cytometry  data  (see,  e.g,  [7,  20,  32]; 
the  problem  is  discussed  at  greater  length  in  [47,  Chs.  1,4]). 

From  Equation  (4.3),  one  sees  that  the  variance  of  the  data  is  not  constant  (as  is  the  case  for  data  described 
by  constant  variance  models),  nor  does  it  scale  with  the  square  of  the  magnitude  of  the  model  (as  is  the  case  for 
data  described  by  constant  coefficient  of  variance  models).  See  [13,  Ch.  3].  Rather,  the  variance  grows  linearly 


I[h](tj,zk;q0) 


(4.3) 


actually  sampled  (as  estimated  by  the 
random  variable 

It  follows  that 

Nl 

=  M 
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with  the  magnitude  of  the  model  solution.  This  relationship  has  been  shown  to  accurately  describe  CFSE  data 
[5,  47],  and  we  can  use  this  relationship  to  establish  a  parameter  estimation  framework. 


4.2  Parameter  Estimation 


The  goal  of  the  parameter  estimation  problem  is  to  find  an  estimate  for  the  parameter  <fo  which  is  assumed  to 
generate  the  data  (neglecting  model  misspecification) .  In  the  above  statistical  model,  we  must  also  estimate  the 
nuisance  parameters  A  =  { Ay } .  Given  a  collection  of  random  variables  Nj  distributed  according  to  (4.3),  one  may 
define  the  estimators 


( qwLS ,  A wls)  =  arg  ^  min  J((q,  A)|{iV^})  =  arg  _  nrin 


(7,A)eQxA 


(«,A)eQxA 


j,k 


(4.4) 


where  the  weights  are  chosen  in  a  manner  that  accounts  for  the  assumed  reliability  of  each  measurement  (see 
below) .  Experimental  data  is  considered  as  a  set  of  realizations  ni  of  the  random  variables  Nj ;  these  are  used  to 
obtain  the  estimates 


(q,  A)  =  arg  min  J((q,  A)|{n^.})  =  arg  min  ^ 


(A  jl[n](tj,zk)q)-n3k^ 


(q,\)eQxA 


(q,\)eQxA 


(4.5) 


j,k 


We  see  that  the  cost  functional  J,  and  thus  the  estimators,  depends  upon  the  data  random  variables  NJk  and 
thus  on  the  statistical  model  (4.3).  In  order  for  the  estimators  qwLS  and  A  wls  to  be  asymptotically  optimal, 
the  weights  wj  must  be  chosen  to  match  the  variance  of  the  random  variables  NJk  [44,  45] , 


im.  = 


Xjf-A n](tj,4;«)),  I[n](tj,z?k\q0)  >  I* 
I[n](tJ,zJk;q0)  <  I* 


\  B_  T* 


(4.6) 


The  cutoff  value  I*  >  0  is  determined  by  the  experimenter  so  that  the  resulting  residuals  appear  random.  In  the 
work  that  follows,  I*  =  200.  The  values  of  B  and  bj  are  known  from  the  experiment.  Notice  that  the  computation 
of  the  weights  (4.6)  depends  upon  the  value  of  the  ‘true’  parameter  qg  as  well  as  on  the  nuisance  parameters  A. 
As  a  result,  one  must  use  an  iterative  estimation  procedure  [13,  19]. 

Traditionally,  it  has  been  assumed  that  A j  =  1  for  all  j  [7,  8,  9,  12,  35]-that  is,  that  there  is  no  scaling  error. 
While  this  assumption  is  obviously  violated  by  some  data  sets  [7,  20,  32,  47]  it  is  not  clear  that  the  incorpo¬ 
ration  or  omission  of  the  nuisance  parameters  will  have  a  significant  effect  on  the  estimation  of  parameters.  In 
practice,  the  nuisance  parameter  vector  must  be  estimated  in  conjunction  with  the  model  parameter  vector  in  a 
two-stage  process.  Unfortunately  two-stage  estimation  may  cause  some  parameters  of  the  mathematical  model 
to  become  unidentifiable  (or,  at  the  very  least,  the  variance  of  the  estimators  for  certain  parameters  may  increase 
dramatically).  For  this  report,  it  will  be  assumed  that  A  j  =  1  for  all  j  and  the  nuisance  parameters  will  not 
be  estimated.  As  will  be  seen  in  Section  5,  the  available  data  is  well-described  by  the  mathematical  model  even 
under  this  simplified  assumption.  We  consider  the  following  estimation  algorithm: 


1.  Set  l  =  0.  Obtain  an  initial  estimate  (the  ordinary  least  squares  estimate)  by  solving  (4.5)  with  A  =  (1, . . . ,  1) 
and  wk  =  1  for  all  j  and  k , 

<7(0)  =  arg mm ^  (j[h](tj ,  zk\  q)  -  n*.)  ; 

j,k 

2.  Compute  the  weights  wJk  for  each  j  and  k  according  to  (4.6)  with  q0  replaced  by  and  with  A  =  (1, . . . ,  1); 

3.  At  iteration  £,  compute  q^  according  to  (4.5)  with  the  current  weights,  and  with  A  =  (1, . . . ,  1); 
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4.  Update  the  weights  again  according  to  (4.6),  now  with  q0  replaced  by  (and  A  =  (1, . . . ,  1)  still);  increment 

t’, 

5.  Repeat  steps  3  and  4  until  convergence  is  obtained. 

5  Results 

Twelve  models  of  cyton  dynamics  are  considered  in  Table  1  and  two  methods  of  estimating  the  statistical  moments 
of  the  autofluorescence  distribution(s)  as  proposed  in  Section  3.4  are  used.  In  order  to  determine  which  of  these 
model  formulations  best  describes  the  available  data,  we  need  mathematical  tools  for  rigorous  model  comparison. 
These  tools  are  explored  in  Section  5.1  below  and  a  best- fit  model  is  selected.  The  fit  of  this  model  to  the  data, 
as  well  as  the  statistical  model  of  the  data,  are  then  discussed.  Finally,  the  dynamic  responsiveness  of  the  cells 
from  the  experimental  data  is  analyzed. 

5.1  Model  Comparison 

From  Section  3.3,  it  is  clear  that  some  of  the  models  of  cyton  dynamics  are  refinements  of  other  models  (in  the 
sense  that  the  more  complex  model  includes  all  possible  solutions  of  the  simpler  model).  For  instance,  the  basic 
cyton  model  (Model  1)  can  be  considered  as  a  refinement  of  a  cyton  model  for  which  it  is  assumed  ipoit)  =  0 
(Model  6).  This  is  because,  for  appropriate  choices  of  the  parameters  E[Tyte]  and  S'TDfT*®],  Model  1  is  exactly 
equivalent  to  Model  6.  In  such  a  case,  it  is  clear  that  the  more  complex  model  must  result  in  a  minimized  cost 
functional  at  least  as  low  as  that  of  the  simpler  model;  yet  the  more  complex  model  has  more  parameters,  and 
one  must  consider  the  possibility  of  overparameterization  against  this  decrease  in  the  least  squares  cost.  To  this 
end,  results  from  the  asymptotic  theory  of  least  squares  estimators  can  be  extended  [5]  to  weighted  least  squares 
estimators  such  as  (4.4). 

Assume  (without  loss  of  generality)  that  the  nuisance  parameters  A  are  known.  Let  qwLS  be  defined,  as 
above,  as  the  minimizer  over  the  admissible  parameter  space  Q  of  the  least  squares  cost  function  J((q. ,  A)|{A^}). 
Now  define  qwLS  to  be  the  minimizer  of  J((q,  A)|{AT^})  over  the  restricted  parameter  set  Qh ,  where  Qh  =  {q  £ 
Q\Hq  =  h}  for  some  linear  function  H  of  rank  r  and  a  vector  h  of  size  r  x  1.  (For  nonlinear  restrictions  on  the 
parameter  space,  a  locally  equivalent  condition  can  be  derived  from  the  first  order  linearization  of  the  nonlinear 
constraint;  see  [5].)  In  such  situations,  the  model  comparison  problem  can  be  recast  as  one  of  hypothesis  testing. 
Consider  the  null  and  alternative  hypotheses, 


Ho  ■  q  G  Qh 
Ha  :  q  Qh- 

Then  under  fairly  general  conditions  (see  [5])  and  assuming  the  null  hypothesis  is  true,  the  test  statistic 

n  ^J((qwLS,  A)|{A^})  -  J{{qwLS,  A)|{A^})) 

lJ yi  —  — > 

J((qwLs,  A)|{JVj£}) 

(where  n  is  the  total  number  of  data  points)  is  asymptotically  a  chi-square  random  variable  with  r  degrees  of 
freedom,  Un  ~  y2(r).  Thus  given  the  data  {n3k}  as  realizations  of  the  random  variables  Nk,  one  obtains  a 
realization  un  of  Un  which  can  be  used  to  assess  the  likelihood  that  the  decrease  in  cost  associated  with  the 
unrestricted  parameter  space  is  the  result  of  chance  (see  [13,  Sec.  3.5]).  The  complete  conditions  under  which  XJn 
is  asymptotically  distributed  as  a  chi-square  random  variable,  as  well  as  a  proof  of  the  result,  can  be  found  in  [5] 
and  the  references  therein. 

For  comparison  among  models  which  are  not  refinements,  one  can  use  information  theoretic  criteria  such  as 
Akaike’s  Information  Criterion  (AIC).  From  (4.3),  it  follows  that  the  scaled  residuals 

=  A j I [n] (tj ,Zk',q)  —  N l  ^  ^ 

p 
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are  independent  and  normally  distributed  with  constant  variance  (for  all  k  and  j).  Then  the  AIC,  which  is  the 
expected  value  of  the  relative  Kullback-Leibler  distance  for  a  given  model  [16],  is 

ku = „  log  +  2p  (5.2, 

where  p  is  the  dimension  of  the  space  Q  for  the  particular  model  of  interest.  Because  Kn  provides  information 
regarding  the  relative  Kullback-Leibler  distance,  AIC  values  have  meaning  only  in  comparison  to  one  another. 
The  model  which  results  in  the  lowest  AIC  value  is  the  most  likely  model  for  the  particular  data  set  being 
investigated.  We  note  that  a  direct  comparison  of  the  costs  (in  Tables  3  and  4)  may  not  be  reasonable  due  to 
the  fact  that  the  AIC  values  also  must  take  into  account  the  number  p  of  parameters  estimated  in  a  given  model 
(see  (5.2)).  A  derivation  of  the  AIC  as  well  as  numerous  examples  can  be  found  in  [16]. 

With  these  tools,  we  are  now  ready  to  compare  the  models  suggested  in  Sections  3.3  and  3.4.  The  minimized 
costs  J{{qwLSi^j))  for  each  donor,  cell  type,  and  model  are  summarized  in  Tables  3  and  4.  Table  3  contains 
the  results  when  the  autofluorescence  distribution  is  assumed  to  be  time-invariant  and  its  statistical  moments  are 
estimated  within  the  least  squares  framework  summarized  in  Section  4;  Table  4  contains  the  results  when  the 
autofluorescence  distribution  is  estimated  at  each  measurement  time  using  the  method  of  moments. 


CD4  T  Cells 

CD8  T  Cells 

Donor  1 

Donor  2 

Donor  1 

Donor  2 

Model  1 

6.0547x  104 

11.812xl04 

7.9383x  104 

20.348  xlO4 

Model  2 

5.4212x  104 

5.1257xl04 

2.8669x  104 

4.5650xl04 

Model  3 

6.0344x  104 

11.812xl04 

8.2439x  104 

20.348  xlO4 

Model  4 

6.1677x  104 

11.863xl04 

7.9383x  104 

20.348  xlO4 

Model  5 

5.4212x  104 

5.1257xl04 

2.7963x  104 

4.5651xl04 

Model  6 

5.4212x  104 

5.1257xl04 

2.8670x  104 

4.5650xl04 

Model  7 

4.3175x  104 

9.4856  xlO4 

6.1827x  104 

18.761xl04 

Model  8 

3.4376x  104 

4.2341  xlO4 

2.8038x  104 

3.9553xl04 

Model  9 

4.4006x  104 

9.5017xl04 

6.0628x  104 

19.563xl04 

Model  10 

4.3196x  104 

9.4931  xlO4 

6.2538x  104 

18.761xl04 

Model  11 

3.4375x  104 

4.2346  xlO4 

2.8004x  104 

3.9551xl04 

Model  12 

3.4376x  104 

4.8931xl04 

2.8038x  104 

3.9554xl04 

Table  3:  Minimized  weighted  least  squares  costs  J((q,  A))  for  various  cyton  model  parameterizations  (see  Table 
1).  Autofluorescence  estimated  as  a  time- invariant  lognormal  density  function  by  least  squares  fit  to  data. 


Of  the  12  models  of  cyton  dynamics  considered,  Model  12  is  generally  selected  as  the  best  model.  The  only 
exception  is  for  Donor  2  CD4  T  cells  when  estimating  a  time-invariant  autofluorescence  distribution  using  least 
squares.  In  this  case,  Model  8  is  narrowly  selected  by  the  AIC  (K  =  9525.77  compared  to  K  =  9525.98),  although 
AIC  differences  less  than  2  are  generally  not  considered  significant  [16].  On  the  other  hand,  the  model  comparison 
test  statistic  (Model  8  is  a  refinement  of  Model  12)  is  un  =  5.7992,  so  that  one  would  reject  the  null  hypothesis 
(Model  12)  only  at  confidences  less  than  87.82%,  which  is  lower  than  typical  thresholds  for  hypothesis  testing. 
Thus  the  results  for  this  particular  data  set  are  ambiguous.  It  should  be  acknowledged  that  Model  8  and  Model 
12  are  quite  similar  (Model  8  is  the  generalization  of  Model  12  allowing  for  undivided  cell  death),  so  that  the 
distinction  between  the  two  is  quite  small.  It  seems  safe  to  consider  Model  12  to  be  the  most  parsimonious  model 
of  the  data  for  both  donors  and  for  both  cell  types. 

A  comparison  of  the  two  methods  of  treating  autofluorescence  is  not  straightforward.  On  one  hand,  the 
minimized  costs  shown  in  Tables  3  and  4  are  directly  comparable  in  the  sense  that  the  costs  do  not  include  any 
measure  of  cost  associated  with  the  autofluorescence  data,  whether  or  not  that  data  is  used.  On  the  other  hand, 
the  estimation  of  a  time-invariant  autofluorescence  distribution  does  not  make  any  use  of  the  autofluorescence 
data.  From  this  perspective  the  two  methods  of  describing  autofluorescence  are  associated  with  distinct  collections 
of  data,  even  if  the  histograms  of  interest  (e.g.,  those  shown  in  Figures  1  and  2)  are  the  same.  It  is  also  not  clear 
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CD4  T  Cells 

CD8  T  Cells 

Donor  1 

Donor  2 

Donor  1 

Donor  2 

Model  1 

7.1255x  104 

14.302xl04 

8.5425x  104 

23.068xl04 

Model  2 

5.8092x  104 

5.4022  xlO4 

3.9701  xlO4 

5.8327xl04 

Model  3 

7.1097x  104 

14.313xl04 

8.4900x  104 

23.503xl04 

Model  4 

7.2613x  104 

14.344xl04 

8.5434x  104 

23.068xl04 

Model  5 

5.8092x  104 

5.4020  xlO4 

3.9709x  104 

5.8309xl04 

Model  6 

5.8092x  104 

5.4041  xlO4 

3.9701  xlO4 

5.8327xl04 

Model  7 

4.3392x  104 

11.098xl04 

6.7689x  104 

20.965xl04 

Model  8 

3.3741  xlO4 

4.4724x  104 

3.7443x  104 

5.2944xl04 

Model  9 

4.3418x  104 

11.098xl04 

6.8680x  104 

20.952xl04 

Model  10 

4.3392x  104 

11.098xl04 

6.7689x  104 

20.965xl04 

Model  11 

3.3741  xlO4 

4.4730xl04 

3.7436x  104 

5.2944xl04 

Model  12 

3.3741  xlO4 

4.4724x  104 

3.7443x  104 

5.2944xl04 

Table  4:  Minimized  weighted  least  squares  costs  J((q,  A))  for  various  cyton  model  parameterizations  (see  Table 
1).  Autofluorescence  estimated  using  the  method  of  moments  at  each  measurement  time. 


whether  the  failure  of  the  minimized  cost  functionals  to  reflect  the  costs  associated  with  the  autofluorescence  data 
is  a  strength  or  a  weakness  of  the  current  approach.  When  such  data  is  available  (as  it  is  here)  it  seems  that  it 
would  make  for  a  useful  comparison.  However,  the  collection  of  such  data  requires  additional  experimental  setup, 
and  this  data  itself  is  only  interesting  to  the  extent  that  it  helps  to  describe  the  dynamics  of  cellular  division  and 
death  as  observed  in  the  histogram  profiles  of  labeled  cells. 

Comparing  the  two  approaches  strictly  in  terms  of  accuracy  in  describing  the  histogram  profiles  of  labeled 
cells  (that  is,  comparing  the  minimized  costs  in  Tables  3  and  4),  the  results  depend  upon  cell  type.  For  CD8  T 
cells,  there  is  a  clear  advantage  in  describing  the  autofluorescence  distribution  as  time-invariant  and  estimating 
the  moments  of  the  distribution  in  a  least  squares  framework.  For  CD4  T  cells,  the  results  are  less  clear.  For 
Donor  1,  there  is  a  very  small  improvement  (among  the  more  accurate  models)  in  using  the  method  of  moments 
to  estimate  the  autofluorescence  distribution.  For  Donor  2,  the  method  of  moments  works  best  for  Model  12,  but 
not  for  Model  8  (which  is  the  AlC-selected  model  when  autofluorescence  is  estimated  by  least  squares).  Again, 
the  difference  is  very  small.  The  analysis  presented  in  the  remainder  of  this  document  is  based  on  results  obtained 
with  Model  12  with  a  time- invariant  autofluorescence  distribution  estimated  in  a  least  squares  framework. 

5.2  Analysis  of  the  Mathematical  and  Statistical  Models 

The  fit  of  the  mathematical  model  to  data  for  both  donors  is  summarized  in  Figures  6  (CD4  T  cells)  and  7  (CD8 
T  cells).  The  figures  show  the  best-fit  model  solution  of  the  mathematical  model  in  comparison  to  the  data;  the 
shaded  region  indicates  the  expected  level  of  ‘noise’  in  the  data  as  a  result  of  the  measurement  process.  That  is, 
if  the  calibrated  mathematical  model  were  perfectly  specified  [E[Nk]  =  I[h](tj,  zy,  <?)),  the  data  would  oscillate 
randomly  around  the  model  solution.  Since  the  data  are  assumed  to  be  normally  distributed  (see  Section  4),  the 
4-standard-deviation  region  highlighted  should  contain  99.9%  of  the  data  points. 

Overall,  the  fit  to  data  is  good.  For  Donor  1  CD4  T  cells,  the  estimated  mean  and  standard  deviation  of 
the  autofluorescence  distribution  are  E[xa]  =  372.42,  STD[xa]  =  220.08.  For  Donor  2  CD4  T  cells  the  estimates 
are  E[xa]  =  543.76,  STD[xa\  =  251.90.  For  CD8  T  cells  the  estimates  are  E[xa]  =  658.35,  STD[xa\  =  319.79 
and  E[xa]  =  542.76,  STD[xa]  =  271.25,  for  Donors  1  and  2,  respectively.  These  numbers  are  within  reason  when 
compared  to  measured  autofluorescence  data  (Figure  5). 

The  primary  shortcoming  of  the  statistical  model  is  the  assumption  of  correct  specification,  i.e. ,  E[Nk\  = 
I[h](tj,  Zk',q )•  There  are  clearly  instances  in  Figures  6  and  7  where  the  data  is  not  centered  around  the  mathemat¬ 
ical  model,  indicating  that  the  mathematical  model  is  systematically  in  error.  As  a  result,  the  shaded  regions  do 
not  contain  99.9%  of  the  data  points.  From  one  perspective,  this  shortcoming  of  the  statistical  model  is  entirely 
mathematical-  if  an  improved  mathematical  model  were  available,  then  it  is  possible  that  the  assumption  of  correct 
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Approximated  Initial  Condition,  t  =23.5  hrs 
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Figure  6:  Calibrated  model  with  CD4  T  cell  data  for  a  particular  set  of  triplicates  from  Donor  1  (left)  and  Donor  2 
(right).  Shaded  regions  indicated  a  4  standard  deviation  confidence  region  computed  according  to  the  theoretical 
statistical  model  (4.3),  assuming  the  model  is  correctly  specified. 
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Approximated  Initial  Condition,  t  =23.5  hrs 
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Figure  7:  Calibrated  model  with  CD8  cell  data  for  a  particular  set  of  triplicates  from  Donor  1  (left)  and  Donor  2 
(right).  Shaded  regions  indicated  a  4  standard  deviation  confidence  region  computed  according  to  the  theoretical 
statistical  model  (4.3),  assuming  the  model  is  correctly  specified. 
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specification  would  be  accurate.  Alternatively,  it  is  possible  to  consider  a  more  general  statistical  model  which 
directly  considers  the  effects  of  misspecification,  for  instance  by  assuming  an  autoregressive  error  structure  [24]. 
This  may  provide  a  more  accurate  measure  for  model  comparison  in  the  absence  of  a  more  accurate  mathematical 
model. 

It  is  also  assumed  in  the  statistical  model  that  V  ar[Nj]  oc  I[h](tj,Zk',q)-  Traditionally,  the  accuracy  of  this 
assumption  is  examined  by  residual  plots  [13,  Ch.  3].  For  instance,  the  modified  residuals  (5.1)  should  be  randomly 
distributed  when  plotted  against  the  magnitude  of  the  model  solution.  However,  this  analysis  is  premised  upon 
the  assumption  of  correct  specification,  and  thus  cannot  be  performed  here.  For  other  data  sets,  it  has  been 
shown  that  the  statistical  model  presented  in  this  document  does  accurately  describe  the  variance  in  CFSE-based 
flow  cytometry  data  [5,  47]. 

In  spite  of  these  shortcomings,  it  should  be  emphasized  that  the  fit  of  the  model  is  quite  good  for  both  donors 
and  cell  types.  As  such,  we  proceed  to  analyze  the  dynamic  responsiveness  of  the  measured  cells  in  the  current 
modeling  framework. 

5.3  Analysis  of  Dynamic  Responsiveness 

From  the  calibrated  mathematical  model,  one  can  compute  the  probability  density  functions  <pi(t)  and  'fi(t)  from 
which  the  times  to  divide  and  die  are  assumed  to  be  drawn  in  the  cyton  model  of  cell  division.  One  can  also 
summarize  the  division  destiny  of  each  population  of  cells.  These  are  summarized  graphically  in  Figure  8  for  CD4 
T  cells  and  Figure  9  for  CD8  T  cells. 

Notice  that  the  cytons  4>o(t)  have  been  truncated  to  the  left  at  t  =  to  =  23.5.  This  is  because  no  information 
is  available  before  the  first  measurement  time;  the  assumption  that  the  density  functions  fa  can  be  described  by 
a  weighted  sum  of  lognormal  densities  does  not  require  that  the  support  of  fa  be  contained  in  the  region  t  >  to, 
so  this  condition  must  be  additionally  imposed  (see  the  Appendix).  This  gives  the  impression  in  Figures  8  and 
9  that  some  fraction  of  cells  will  begin  to  divide  immediately  following  the  first  measurement  time.  Though  this 
may  be  true,  it  is  beyond  the  capabilities  of  the  current  modeling  framework  to  determine.  This  is  because  the 
estimated  parameters  are  only  unique  up  to  the  numbers  of  cells  predicted  to  divide  and  die  between  measurement 
times.  In  other  words,  if  {tj}  is  a  collection  of  measurement  times,  the  model  provides  meaningful  estimates  of 
the  quantities 


the  fraction  of  cells  from  the  initial  population  estimated  to  have  completed  the  first  division  between  measurement 
times  tj  and  tj+ These  quantities  (converted  to  percentages)  are  superimposed  on  the  graphs  of  the  curves 
Fofa(t)  in  Figures  8  and  9.  Thus,  although  the  current  modeling  framework  does  not  provide  information 
regarding  when  cells  will  begin  to  divide  it  does  provide  meaningful  information  on  the  distribution  of  times  to 
first  division.  Unsurprisingly,  this  information  is  constrained  by  the  frequency  at  which  the  population  is  measured. 

We  see  that  CD4  T  cells  from  Donor  1  reach  their  first  division  more  rapidly  and  in  greater  quantity  than 
CD4  T  cells  from  Donor  2.  By  the  end  of  the  experiment,  75.34%  of  the  initial  population  of  CD4  T  cells  from 
Donor  1  had  divided  in  response  to  stimulus,  compared  to  69.50%  for  Donor  2.  CD8  T  cells  from  Donor  1  likewise 
respond  more  rapidly  and  in  greater  quantity  than  those  from  Donor  2.  By  the  end  of  the  experiment,  88.44%  of 
the  initial  population  of  cells  from  Donor  1  had  divided,  while  only  63.94%  of  those  from  Donor  2  had  divided. 
Comparing  CD4  T  cells  and  CD8  T  cells  within  the  same  donor,  we  observe  that  CD8  T  cells  complete  their  first 
division  more  quickly  than  CD4  T  cells. 

For  cells  having  already  completed  at  least  one  division,  CD4  T  cells  from  Donor  1  are  estimated  to  divide 
more  slowly  and  with  greater  variation  across  the  population  (E[Tfvv]  =  21.2,  STD[Tfm]  =  19.5)  compared 
with  those  from  Donor  2  (. E[Tflv ]  =  11.3,  STD[Tfzv]  =  5.7).  Similar  behavior  is  observed  for  CD8  T  cells 
(. E[T?iv ]  =  11.2,  STD[Tfiv]  =  4.6  for  Donor  1;  E[Tfiv }  =  9.4,  STD[Tfiv]  =  3.0  for  Donor  2). 

The  analysis  of  cell  death  is  complicated  by  several  factors.  While  the  cyton  density  functions  fa  and  fa 
are  assumed  to  be  identical  for  i  >  1  (that  is,  for  all  cells  having  divided  at  least  once),  the  fraction  F,;  of 
progressing  cells  varies  from  one  generation  to  the  next  according  to  Equation  (3.12).  The  fraction  of  cells  which 
are  non-progressors  dies  according  to  the  density  function  fa,  but  will  not  divide.  Thus  the  division  destiny  for 
the  population  of  cells  must  be  taken  into  account  when  interpreting  cell  death.  Simultaneously,  even  among  cells 
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Undivided  Cells:  Time  to  First  Division  Undivided  Cells:  Time  to  First  Division 


Division  Destiny  Division  Destiny 


Divisions  Divisions 


Figure  8:  Comparison  of  estimated  CD4  T  cell  division  dynamics  between  Donor  1  (left)  and  Donor  2  (right). 
Top:  Probability  density  function  <po(t)  for  time  to  first  division  for  initially  undivided  CD4  T  cells,  scaled  by 
the  initial  progressor  fraction  Fq  of  activated  cells.  Percentages  indicate  the  fraction  of  undivided  cells  which  will 
have  entered  their  first  division  by  the  next  measurement  time  (vertical  dashed  lines).  On  average,  cells  for  Donor 
1  complete  their  first  division  more  rapidly  than  those  for  Donor  2;  cells  from  Donor  1  are  also  more  likely  to  have 
divided  in  response  to  stimulus  before  the  end  of  the  experiment.  CD4  T  cells  from  both  donors  are  estimated 
to  respond  more  slowly  and  in  greater  frequency  than  CD8  T  cells  (compare  Figure  9).  However  the  total  CD  8 
T  cell  response  for  Donor  1  is  greater  than  the  CD4  T  cell  response  while  the  total  CD4  and  CD8  responses  are 
comparable  for  Donor  2.  Middle:  Probability  density  functions  for  time  to  subsequent  division  or  time  to  die 
(inverted)  for  CD4  T  cells  having  completed  at  least  one  division.  Cells  from  Donor  2  divide  slightly  more  rapidly 
and  more  synchronously  than  those  from  Donor  1.  Bottom:  Division  destiny,  indicating  the  average  number 
of  divisions  undergone  by  cells  initially  in  the  population  (at  t  =  to).  The  fraction  of  cells  with  division  destiny 
equal  to  zero  estimates  the  relative  abundance  of  cells  which  will  not  become  activated  to  divide. 
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Undivided  Cells:  Time  to  First  Division  Undivided  Cells:  Time  to  First  Division 


Figure  9:  Comparison  of  estimated  CD8  T  cell  division  dynamics  between  Donor  1  (left)  and  Donor  2  (right). 
Top:  Probability  density  function  <f>o(t)  for  time  to  first  division  for  initially  undivided  CD8  T  cells,  scaled  by 
the  initial  progressor  fraction  Fq  of  activated  cells.  Percentages  indicate  the  fraction  of  undivided  cells  which  will 
have  entered  their  first  division  by  the  next  measurement  time  (vertical  dashed  lines).  On  average,  cells  for  Donor 
1  complete  their  first  division  more  rapidly  than  those  for  Donor  2  and  cells  from  Donor  1  are  more  likely  to  have 
divided  in  response  to  stimulus  before  the  end  of  the  experiment.  Middle:  Probability  density  functions  for  time 
to  subsequent  division  or  time  to  die  (inverted)  for  CD8  T  cells  having  completed  at  least  one  division.  Cells 
from  Donor  2  divide  slightly  more  rapidly  and  more  synchronously  than  those  from  Donor  1.  Bottom:  Division 
destiny,  indicating  the  average  number  of  divisions  undergone  by  cells  initially  in  the  population  (at  t  =  to)-  The 
fraction  of  cells  with  division  destiny  equal  to  zero  estimates  the  relative  abundance  of  cells  which  will  not  become 
activated  to  divide. 


27 


which  would  be  progressors,  cell  death  is  a  possibility  if  that  cell’s  time-to-die  (sampled  according  to  the  density 
ipi  is  less  than  time-to-divide  (sampled  according  to  the  density  (pi). 

Based  upon  the  densities  (pi  and  ipi  ( i  >  1)  in  Figures  8  and  9  and  ignoring  division  destiny,  it  would  appear 
that  cells  from  Donor  1  (both  CD4  T  and  CD8  T)  are  more  likely  to  die  relative  to  cells  taken  from  Donor  2.  Yet 
when  the  numbers  of  dying  cells  are  computed  (Figure  10),  we  find  that  this  hypothesis  holds  for  CD4  T  cells 
(cells  from  Donor  1  are  slightly  more  prone  to  die)  but  completely  fails  for  CD8  T  cells.  When  division  destiny  is 
taken  into  account,  we  see  that  CD4  T  cells  from  Donor  1  are  estimated  to  have  a  much  narrower  division  destiny 
than  CD4  T  cells  from  Donor  2.  As  a  result,  fewer  cells  progress  to  high  division  number  (i  >  6),  and  these 
non-progressors  begin  to  die.  The  result  is  that  more  cells  are  observed  with  high  division  number  for  Donor  2, 
and  the  expansion  of  the  population  of  cells  over  the  course  of  the  experiment  is  greater  for  Donor  2.  For  CD8  T 
cells,  the  estimated  division  destiny  for  Donor  2  is  narrow  but  has  a  high  mean.  Coupled  with  the  faster  rate  of 
division  (middle-right  panel,  Figure  9),  we  see  that  most  CD8  T  cells  from  Donor  2  proceed  through  large  number 
of  divisions,  after  which  the  population  has  reached  its  division  destiny  and  the  size  of  the  population  reaches  a 
maximum  before  cells  begin  to  die. 

Thus,  in  effect,  division  destiny  imposes  a  limit  on  the  degree  to  which  a  population  of  cells  can  expand  in 
response  to  stimulus.  For  CD4  T  cells  from  Donor  1,  the  population  appears  to  have  reached  its  maximum 
expansion  just  as  the  experiment  ends,  expanding  by  a  total  factor  of  6.01  (ratio  of  maximum  population  size  to 
initial  size).  For  CD4  T  cells  from  Donor  2,  the  expansion  factor  at  the  end  of  the  experiment  is  9.75  and  still 
rising.  For  CD8  T  cells,  Donor  1  expands  by  a  factor  of  20.88  (and  rising)  by  the  end  of  the  experiment,  while 
for  Donor  2  the  population  expanded  by  a  factor  of  23.55  before  it  began  to  contract. 

Therefore  we  see  that  the  cytons  </>,;  and  i pi  provide  meaningful  information  regarding  the  times  at  which  cells 
will  divide  or  die,  but  this  information  must  be  carefully  interpreted  with  respect  to  division  destiny.  This  can  be 
accomplished  by  reconstructing  the  population  generation  structures  for  viable  and  dead  cells  (as  in  Figure  10). 
Then,  one  can  make  deductions  concerning  the  viability  of  the  populations  of  cells  by  analyzing  the  numbers  of 
cells  as  described  above. 
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Figure  10:  Cell  numbers  and  population  generation  structure  for  CD4  T  cells  (top)  and  CD8  T  cells  (bottom) 
for  Donor  1  (left)  and  Donor  2  (right).  Shaded  areas  represent  model  generated  numbers  while  dots  represent 
experimental  data  values.  Total  numbers  of  dead  cells  (with  generation  structure)  are  shown  inverted.  Note 
that  the  numbers  of  dead  cells  are  cumulative  and  does  not  reflect  any  decrease  in  numbers  of  dead  cells  which 
would  be  associated  with  disintegration.  While  it  is  difficult  to  determine  the  numbers  of  dying  cells  from  cyton 
graphics  alone  (Figures  8  and  9),  one  can  clearly  compare  the  relative  frequency  of  cell  death  between  donors 
and  cell  types  using  these  reconstructions  of  the  population  behavior.  For  both  donors  and  cell  types,  cell  death 
is  estimated  to  be  negligible  until  3-4  days  after  stimulation  (not  considering  any  cells  which  die  before  the  first 
measurement  is  taken).  After  this  time,  most  cells  appear  to  have  reached  their  division  destinies  (see  Figures  8 
and  9,  bottom)  after  which  time  they  begin  to  die. 
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6  Concluding  Remarks 

In  this  document,  we  have  described  and  analyzed  a  recent  class  of  mathematical  models  which  combines  the 
cyton  model  of  population  generation  structure  with  a  mass-conservation  model  of  label  dynamics.  Unlike  previous 
label-structured  models,  the  new  class  of  models  describes  the  processes  of  cellular  division  and  death  in  intuitive 
terms  which  are  relatable  to  important  biological  features.  Significantly,  because  the  new  models  can  be  fit 
directly  to  CFSE  histogram  data,  it  is  possible  to  consider  the  statistical  properties  of  such  data.  From  these 
properties  and  under  mild  assumptions,  a  statistical  model  of  the  data  has  been  derived  and  incorporated  into  a 
least  squares  parameter  estimation  framework.  Using  this  framework,  various  models  selected  from  the  new  class 
of  models  were  fit  to  experimental  data  and  compared.  The  best-fitting  model  has  been  observed  to  accurately 
describe  the  behavior  of  both  CD4  T  cells  and  CD8  T  cells  acquired  from  two  healthy  donors  and  stimulated  to 
divide  with  PHA. 

Of  those  models  tested,  the  selected  model  (Model  12)  features  a  bimodal  distribution  of  times  to  first 
division  and  ignores  cell  death  for  undivided  cells.  From  the  distribution  of  times  to  first  division,  it  is  possible 
to  compute  the  fraction  of  cells  completing  the  first  division  between  each  set  of  measurements.  Distributions 
for  cell  division  and  death  for  cells  having  already  completed  at  least  one  division  are  assumed  to  be  described 
by  lognormal  density  functions.  Though  the  best-fit  mathematical  model  is  observed  to  be  accurate,  there  is 
some  room  for  improvement.  To  this  end,  the  modeling  framework  presented  here  is  readily  generalizable;  any 
distributions  of  times  to  divide  and/or  die  which  admit  density  functions  can  be  tested.  Moreover,  because  the 
mathematical  solution  is  separable  (see  Proposition  3.2)  the  cyton  model  of  population  generation  structure  can 
be  replaced  by  any  other  model  of  cellular  dynamics  with  a  sufficiently  similar  form  (e.g.,  branching  process 
models  [22,  32,  38,  39]).  It  is  possible  that  the  model  may  be  improved  further  by  a  more  detailed  consideration 
of  the  effects  of  cellular  autofluorescence  and  the  changes  of  the  autofluorescence  distribution  over  time. 

The  primary  shortcoming  of  the  statistical  model  is  the  assumption  that  the  model  is  correctly  specified; 
otherwise,  the  statistical  model  has  been  previously  shown  to  correctly  account  for  the  variance  in  histogram 
data  [5,  47].  Thus,  for  a  sufficiently  accurate  mathematical  model,  the  statistical  model  of  CFSE  data  presented 
here  can  be  incorporated  into  a  model  comparison  framework  so  that  alternative  mathematical  descriptions  of 
cell  division  can  be  tested  in  a  manner  that  is  statistically  rigorous.  However,  the  lack  of  inclusion  of  model 
misspecification  in  the  statistical  model  suggests  that  to  use  this  statistical  model  for  computation  of  confidence 
intervals  via  either  asymptotic  theory  or  bootstrapping  may  not  be  appropriate.  Moreover,  the  statistical  model 
given  by  (4.3)  is 

B  i 

Nfc  =  ,  Zk)  +  (A j ,  Zk))* £kj i 

bj 

where  £kj  ~  jV(0,  1),  was  derived  by  considering  a  single  histogram  bin  at  a  single  time.  The  statistical  model  re¬ 
sults  from  repeating  this  derivation  for  each  histogram  bin  at  each  measurement  time.  However,  in  this  derivation, 
the  dependence  of  the  cell  counts  (and  thus  the  probabilities)  on  additional  factors  has  been  completely  ignored. 
The  model  values  I[n](tj,Zk)  will  depend  strongly  on  the  set  of  bins  [zk,  Zk+ i)  used  for  the  histogram  data.  It  can 
be  readily  seen  that  the  level  of  noise  in  the  data  (relative  to  the  magnitude  of  the  data)  increases  as  the  number 
of  bins  increases.  Conversely,  as  fewer  bins  are  used,  the  data  is  effectively  ‘smoothed  out’  or  averaged  and  some 
smaller  features  of  the  population  data  may  be  lost.  Thus  there  seems  to  be  some  optimal  number  of  bins  to 
use  to  represent  the  histogram  data.  On  one  hand,  this  possibility  could  be  assessed  by  trial  and  error  on  the 
number  of  histogram  bins.  Alternatively,  the  statistical  model  might  be  analyzed  and/or  generalized  to  explicitly 
incorporate  the  dependence  of  the  statistical  model  on  the  number  of  histogram  bins,  and  more  importantly,  the 
dependence  of  parameter  confidence  intervals  on  the  number  of  histogram  bins.  Finally,  the  statistical  model 
makes  the  simplifying  assumption  that  the  numbers  of  cells  counted  into  each  distinct  histogram  bin  represents 
an  independent  (from  the  other  bins)  process.  But  this  is  not  true,  for  if  Sj  cells  are  measured  at  time  tj,  then  we 
must  have  the  identity  ]>/,.  Ml  =  Sj.  Thus  the  random  variables  representing  the  numbers  of  cells  Nl  =  S-Mj 

in  the  total  population  counted  into  distinct  bins  are  not  independent.  So,  for  various  reasons  we  cannot  expect 
to  be  able  to  compute  standard  errors  or  confidence  bounds  for  the  estimated  parameters  in  an  unbiased  manner 
[2,  13]. 

The  work  presented  here  demonstrates  how  the  current  modeling  framework  can  be  used  as  a  basis  for 
comparison  between  multiple  donors  and/or  cell  types.  There  are  two  primary  limitations  for  such  comparisons. 
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First,  while  a  comparison  among  donors  and/or  cell  types  may  focus  on  the  differences  in  estimated  moments 
and/or  cyton  density  functions  (such  as  those  shown  in  Figures  8  and  9),  the  information  contained  within  these 
estimated  distributions  is  limited  to  the  chosen  modeling  framework.  Thus,  for  instance,  one  cannot  determine 
the  time  at  which  cells  first  begin  to  divide  only  from  this  information.  Of  course,  more  complex  models  can 
be  incorporated  into  the  current  modeling  framework  if  knowledge  of  this  information  should  prove  necessary. 
Second,  there  has  not  been  (to  our  knowledge)  a  comprehensive  study  of  the  biological  and  experimental  variability 
inherent  in  the  measurement  process.  In  other  words,  it  is  not  known  how  the  behavior  of  cells  from  a  single  donor 
may  vary  from  day  to  day  (if  multiple  blood  samples  are  acquired)  or  from  sample  to  sample  (even  if  acquired  at 
the  same  time). 

Because  the  Malthusian  cell  proliferation  and  death  rates  of  [7,  8]  are  not  necessarily  compatible  with  a 
requirement  of  minimum  cell  cycle  times,  we  have  here  incorporated  and  used  the  cyton  models  of  Section  3.3.  As 
noted  above  (see  (3.10))  this  new  formulation  is  compatible  with  time  dependent  Malthusian  death  rates  ($%{£)  (and 
in  some  cases  with  time  dependent  Malthusian  proliferation  rates  ai(t)).  However,  several  other  generalizations 
of  the  proliferation  and  death  rate  terms  are  immediately  available. 

One  might  consider,  for  example,  the  addition  of  a  second  structure  variable  (say,  volume  or  physiological 
age  [15]),  which  could  be  used  to  enforce  a  minimum  cell  cycle  time  by  requiring  that  cells  progress  from  some 
size  V  to  2V  before  dividing,  at  which  point  two  cells  of  size  V  are  produced.  However,  in  the  absence  of 
additional  observations,  it  is  unclear  what  parameters  (e.g.,  average  rate  of  growth,  or  the  structure  variable  V) 
could  be  estimated  from  CFSE  histogram  data.  Video  microscopy  measurements  by  Hawkins  et  al.  [30]  indicate 
that  average  cell  size  may  be  division  dependent,  and  this  may  complicate  the  inclusion  of  volume  structure. 
Biologically,  it  is  expected  that  apoptosis  occurs  only  at  particular  checkpoints  in  the  cell  cycle  (particularly  if 
external  ‘kill  signals’  are  absent),  so  that  a  generalization  to  volume  structure  (or  any  other  surrogate  for  cell 
cycle  position  or  physiological  age  [15])  may  permit  a  more  accurate  description  of  cell  death.  Still,  it  is  unclear 
what  information  might  be  available  when  considering  only  CFSE  histogram  data.  It  is  possible  that  the  forward 
scatter  (FSC)  of  laser  light  might  be  used  as  some  sort  of  observable  surrogate  for  cell  size,  but  additional  work 
will  be  necessary  to  investigate  this  hypothesis.  However,  a  more  promising  approach  may  involve  use  of  recently 
developed  fluorescence  microscopy  data  (such  as  the  fluorescent  ubiquitination-based  cell  cycle  indicator  (FUCCI) 
in  [15])  to  estimate  probability  density  functions  representing  durations  of  cell  cycle  phases. 

Ideally,  cell  cycle  parameters  (as  represented  in  the  cyton  model)  can  be  related  back  to  more  physi¬ 
cally/experimentally  meaningful  parameters  such  as  the  type  and  strength  of  stimulation,  which  may,  in  turn,  re¬ 
quire  the  translation  of  certain  molecular  pathways  within  individual  cells  into  mathematical  equations/expressions. 
Recent  work  has  indicated  that  the  mechanisms  responsible  for  cell  proliferation  and  death  may  be  mutually  de¬ 
pendent  upon  a  common  molecular  pathway  [21,  46].  As  more  data  becomes  available,  we  hope  to  examine  how 
the  estimated  parameters  change  under  various  experimental  conditions,  with  an  eye  toward  additional  constitu¬ 
tive  relationships  linking  molecular  and/or  subcellular  functions  to  population  dynamics  [17].  In  this  context,  it 
seems  necessary  to  consider  the  extent  to  which  these  functions  and/or  pathways  are  inherited.  Evidence  sug¬ 
gests  that  closely  related  cells  exhibit  strong  correlation  in  times  to  divide  and  some  correlation  in  times  to  die, 
and  that  this  correlation  tends  to  decrease  with  the  number  of  divisions  undergone  [30].  Cells  with  a  common 
precursor  may  also  share  a  common  division  destiny  [30],  which  can  be  altered  by  stimulation  conditions  [48]. 
While  computed  cell  numbers  are  relatively  unaffected  provided  correlation  is  limited  to  cells  having  undergone 
the  same  number  of  divisions  [22,  30,  32],  correlation  between  subsequent  division  of  cells  can  alter  the  dynamics 
predicted  by  a  mathematical  model  [50].  For  large  populations,  this  effect  seems  negligible,  but  may  play  an 
important  role  in  vivo  where  only  a  small  number  of  responding  cells  can  trigger  an  immune  response  [50].  As 
noted  above,  branching  process  models  have  been  formulated  to  account  for  various  levels  of  correlation,  and 
these  models  may  be  incorporated  into  the  compartmental  model  framework  as  described  above. 

In  spite  of  the  limitations  discussed  above,  the  proposed  mathematical  and  statistical  framework  represents 
a  positive  step  toward  a  more  comprehensive  model  of  cellular  division  as  measured  by  flow  cytometry.  The 
flexibility  of  the  class  of  mathematical  models  combined  with  the  probabilistic  treatment  of  the  data  collection 
process  allows  for  a  rigorous  comparison  between  competing  descriptions  of  cellular  behavior.  As  such,  this 
framework  can  help  to  test  biological  hypotheses  and  serve  as  a  bridge  between  cellular-level  events  and  population- 
level  observations.  This  framework  can  also  be  used  to  study  the  optimal  design  of  experiments;  thus  it  may  be 
possible  to  identify  measurement  times  which  minimize  the  size  of  the  blood  sample  required  while  maximizing  the 
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information  one  can  obtain.  In  the  future,  it  will  be  possible  to  analyze  cellular  behavior  for  donors  in  a  variety 
of  clinical  states  and  thus  to  develop  a  more  complete  understanding  of  infectious  disease,  immunosuppressive 
drug  actions,  etc. 
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A  Comments  on  Numerical  Methods  and  Code 


The  computational  algorithm  for  the  model  (3.9)  is  an  extension  (accounting  for  the  incorporation  of  cyton  dynam¬ 
ics)  of  the  algorithm  originally  proposed  by  Allgower  et  ah,  [27].  Because  the  solution  is  factorable  (Proposition 
3.2),  it  is  possible  to  compute  the  population  generation  structure  ( Ni(t ),  for  i  >  0  and  t  >  0)  independently, 
and  then  to  use  this  information  to  compute  the  population  label  structure  (fb(t,  x),  for  i  >  0,  t  >  0,  and  x  >  0). 
Naively,  one  could  compute  the  functions  Ni(t)  numerically  and  the  functions  fii(t,x)  either  numerically  or  ex¬ 
actly  (depending  on  the  form  of  the  function  v(t)).  One  then  obtains  the  densities  n,;(t,  x)  by  Proposition  3.2  and 
the  densities  hi(t,x)  by  the  convolution  integral  (3.5).  However,  this  naive  approach  is  computationally  intensive 
as  a  result  of  the  convolution.  As  shown  in  [27],  there  is  a  more  efficient  method.  We  first  discuss  the  overall 
computational  scheme  for  the  construction  of  the  population  label  structure,  followed  by  a  detailed  algorithm  for 
the  computation  of  the  population  generation  structure. 


A.l  Computation  of  Population  Label  Structure 

Assume  one  has  already  computed  the  functions  _/V,;(f).  The  solutions  fb(t,  x)  of  (3.4)  can  be  obtained  using  the 
method  of  characteristics, 

hi(t,  x)  =  -A-T  ^2lxe(^o  v(s)ds)'j  .exp  v(s)ds ^  .  (A.l) 


Now,  assume  (for  the  moment)  that  the  initial  label  density  in  the  population  is  lognormally  distributed  with 
parameters  /z o  and  cro  so  that 
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for  x  >  0.  Inserting  this  definition  into  (A.l), 
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where 


Hi(t)  =  —  zlog2  —  /  v(s)ds  +  hq. 
Jt0 


In  other  words,  if  the  initial  label  density  is  lognormally  distributed,  then  the  distribution  of  CFSE  (that  is, 
the  distribution  of  fluorescence  intensity  resulting  from  CFSE)  will  be  lognormally  distributed  at  all  times,  with 
parameters  /z,  (t)  and  <To- 

Now,  assume  more  generally  that  the  initial  condition  can  be  written  as  a  convex  combination  of  lognormal 
density  functions, 

K 

$(a?)  /fclogn(x;  nk,  (crfc)2), 

fc= l 


This  assumption  is  not  overly  restrictive  and  the  initial  condition  (see  the  measurements  collected  on  Day  1  in 
Figures  1  and  2)  can  be  well-approximated  by  such  a  series.  Then  by  the  principle  of  superposition,  Proposition 
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3.2,  and  Equations  (4.1)  and  (A. 2) 


K 


n(t,x)  =  Y^Ni(t)'^2fklogn(x;fi%{t),  (ak)2), 

2=0  k= 1 

where  /if  (f)  is  given  as  above.  To  account  for  the  contributions  of  cellular  autofluorescence,  we  have  from  (3.5) 
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By  assumption,  p(t,  £)  is  itself  a  lognormal  density  function  and  the  integral  above  (for  each  pair  of  values  (*,  k))  is 
the  convolution  of  two  lognormal  density  functions,  which  can  be  accurately  approximated  by  a  lognormal  density 
function  having  a  mean  and  variance  which  is  the  sum  of  the  means  and  variances  of  the  two  density  functions 
in  the  convolution  (see  [27]  and  the  references  therein).  In  other  words 

00  K 

n{t,x)  &YNS)  /fclogn(x;  /tf  (t),  (of  (t))2)  (A.3) 

2  =  0  k=  1 


where 


flk(t)  =  log  (Ek(t)) 


(  STDk(t) 

v  Em 


&H0 

Ejm 

STEPS) 


(STDS) 
V  Ek(t ) 


exp  ^(t)  +  +  E[xa] 

((gl^)--1)  .  exp  (2 rf(t)  +  {crk)2) 


Thus,  given  the  population  generation  structure  iV;(i),  the  values  {/&},  {pk},  and  {(<Jk)2}  which  represent  the 
initial  condition  ^(a;),  and  the  parameters  E[xa]  and  STD[xa]  describing  the  distribution  of  autofluorescence,  one 
can  very  quickly  construct  the  solutions  n(t,  x)  using  (A.3).  Significantly,  this  approximation  to  the  solution  does 
not  involve  any  discretization  in  the  structure  variable  (x  or  x)  so  that  solutions  can  be  evaluated  cheaply  even 
on  a  very  fine  mesh  in  the  structure  variable.  We  find,  in  agreement  with  [27],  that  this  method  of  approximation 
increases  computational  speed  by  several  orders  of  magnitude  over  other  methods  of  solution  (e.g.,  [12]). 

Once  the  label-structured  densities  h,:(t,  x)  have  been  obtained,  we  must  compute  the  cell  counts  I[n](tj,  Zk',  q) 
according  to  (4.2).  The  values  n(t,  x)  can  be  computed  very  cheaply  and  efficiently  by  (A.3)  so  that  a  large  number 
of  evaluations  of  h(t,x)  can  be  used  in  the  approximation  of  the  integral  operator  I[h](tj,  Zk\  q)  with  no  adverse 
effect  on  computational  time.  For  the  results  presented  in  this  manuscript,  the  values  I[n\(tj,Zk',q)  have  been 
approximated  using  two  point  Gauss-Legendre  quadrature  on  each  interval  \zk,  Zk+ 1]- 


A. 2  Computation  of  Population  Generation  Structure 

We  now  discuss  a  computational  scheme  for  a  general  cyton  model,  given  by  Equations  (3.6)  -  (3.8).  Assume 
and  ipi(t)  are  known  functions  of  time  for  all  i  >  0.  Given  initial  and  final  measurement  times  to  and  </,  as 
well  as  a  time  step  size  h,  define  the  number  of  time  steps 

N=t-l 

h 
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and  the  time  grid  points 


tj  =  t0  +  (j  -  l)h,  j  =  1,...,(N+1). 

For  each  j,  the  values  4>o(tj)  and  ipoitj)  can  be  precomputed  for  each  i  >  0  and  stored  in  vectors  of  size  (TV  +  1). 
Similarly,  the  values  ft3  4>o (s)ds  and  ft3  il>o (s)ds  can  be  precomputed  and  stored  in  vectors.  The  integration  can  be 
efficiently  carried  out  using  two-point  Gauss-Legendre  quadrature  on  each  subinterval  of  size  h.  We  have  generally 
found  h  =  1  to  be  a  sufficiently  small  time  step.  The  results  in  this  document  were  all  obtained  with  h  =  0.25. 
Because  precomputation  is  cheap,  computation  of  the  terms  ft3  (j>o(s)ds  and  ft °  ipo( s)ds  using  a  higher-order  rule 
(e.g.,  Gauss-Legendre  quadrature)  allows  for  a  larger  value  of  h  than  would  otherwise  be  acceptable. 

When  the  functions  <f>o(t)  and  ipo{t)  are  parametric  density  functions  (as  is  the  case  in  this  document),  it  is 
possible  that  a  portion  of  the  support  of  the  functions  lies  in  the  half  line  t  <  to.  In  order  for  the  cyton  model  to 
function  properly,  it  must  be  true  that  ft  </>i(s)ds  =  ft  ^i(s)ds  =  1  for  all  i  >  0.  Since  one  does  not  have  any 
information  about  the  behavior  of  the  population  of  cells  prior  to  t  =  to,  the  simplest  method  of  resolving  this 
problem  is  to  truncate  (on  the  left)  the  functions  <f>o(t)  and  if>o(t)  at  t  =  to-  Thus  we  can  define 


4>o(t) 

4>o(t) 


4>o(t) 

1  -  Jo°  0o(s)ds 

V’o  (t) 

1  -  fo°  V’ o(s)ds 


If  one  of  the  denominators  above  is  zero,  we  define  (po(t)  (or  ipo(t))  to  be  identically  zero.  We  will  simply  refer 
to  4>o(t)  and  ipo(t)  without  tildes,  although  it  should  be  understood  that  the  functions  have  been  appropriately 
scaled. 

Given  the  precomputed  vectors  above,  one  can  compute  the  quantities  n^v (tj )  and  n*e(G)  according  to  (3.7) 
for  each  j.  These  equations  can  be  computed  for  all  values  of  j  simultaneously  using  a  single  element-wise  vector 
multiplication.  From  the  quantities  n(jtv(tj)  and  n*e(G)>  one  can  obtain  N0(tj )  using  the  trapezoidal  rule. 

Next,  we  can  define  an  additional  vector  of  time  values,  Sj  =  tj  for  all  j  =  1, . . . ,  ( N  +  1).  Then  the  values  of 
4>i{tj  —  Sfc),  ipi(tj  —  Sk),  fy3  Sk  <f>i(£)dt;  and  f^3  ”k  </>i(£)d£  can  be  precomputed  for  i  >  1  and  stored  in  an  array 

of  size  (N  +  1)  x  (N  +  1).  To  do  so  efficiently,  one  can  simply  compute  /0*J  <pi(s)ds  (and  similarly  for  ipi)  where 
tj  =  tj  —  to  for  each  j  and  store  these  quantities  in  a  vector;  the  (j,  k)  entry  of  each  array  is  the  ( j  —  k  +  1)  entry 
of  the  corresponding  vector,  or  is  zero  if  j  —  k  +  1  <  1.  Note  that  in  this  document,  </>,;(f)  and  4>i(t)  are  identical 
for  all  i  >  1,  so  that  only  4  arrays  are  required  to  store  these  precomputed  values. 

Given  these  precomputed  vectors,  one  can  compute  (recursively  on  i)  the  values  nftv(tj)  and  nfle(tj)  according 
to  (3.8)  for  each  j.  Again,  these  values  can  be  computed  for  all  values  of  j  simultaneously  by  carefully  vectorizing 
the  resulting  operations.  The  terms  n*1)  in  (3.8)  are  vectors  of  size  ( N  +  1),  which  can  be  replaced  by  an 
( N  +  1)  x  ( N  +  1)  array  where  each  row  is  the  vector  of  values  nf^itj).  Then  the  equations  (3.8)  can  be 
computed  using  element-wise  matrix  multiplication,  followed  by  quadrature  (using  the  trapezoidal  rule)  over 
values  of  Sk ■  From  the  quantities  nflv(tj)  and  nfle(tj),  one  can  obtain  Ni(tj)  using  the  trapezoidal  rule. 


A. 3  Inverse  Problem/Parameter  Estimation 

Given  the  forward  solution  constructed  as  described  in  the  previous  subsections,  this  information  must  be  incor¬ 
porated  into  a  computational  scheme  for  the  optimization  problem  (4.5).  This  optimization  is  carried  out  using 
the  BFGS  algorithm  as  implemented  in  the  Matlab  routine  fmincon.  Computations  were  carried  out  on  a  Dell 
Optiplex  990,  running  an  Intel  Core  i7-2600  (4x3.4GHz)  with  2x4BG  RAM  (1333MHz).  The  inverse  problem 
took  an  average  of  6.41  minutes. 
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